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

changed to new state, please report bugs to Luv or Martin

This commit is contained in:
Martin Diehl 2014-07-02 12:27:39 +00:00
parent a68fe7c77b
commit 8fa2dcffbd
21 changed files with 1197 additions and 3893 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

56
code/CPFEM.f90
View File

@ -75,12 +75,10 @@ subroutine CPFEM_initAll(temperature,el,ip)
use FEZoo, only: & use FEZoo, only: &
FEZoo_init FEZoo_init
#endif #endif
#ifdef NEWSTATE
use constitutive_thermal, only: & use constitutive_thermal, only: &
constitutive_thermal_init constitutive_thermal_init
use constitutive_damage, only: & use constitutive_damage, only: &
constitutive_damage_init constitutive_damage_init
#endif
implicit none implicit none
integer(pInt), intent(in) :: el, & ! FE el number integer(pInt), intent(in) :: el, & ! FE el number
@ -108,10 +106,8 @@ subroutine CPFEM_initAll(temperature,el,ip)
call lattice_init call lattice_init
call material_init call material_init
call constitutive_init call constitutive_init
#ifdef NEWSTATE
call constitutive_thermal_init call constitutive_thermal_init
call constitutive_damage_init call constitutive_damage_init
#endif
call crystallite_init(temperature) ! (have to) use temperature of first ip for whole model call crystallite_init(temperature) ! (have to) use temperature of first ip for whole model
call homogenization_init call homogenization_init
call CPFEM_init call CPFEM_init
@ -156,10 +152,6 @@ subroutine CPFEM_init
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase material_phase
#ifndef NEWSTATE
use constitutive, only: &
constitutive_state0
#endif
use crystallite, only: & use crystallite, only: &
crystallite_F0, & crystallite_F0, &
crystallite_Fp0, & crystallite_Fp0, &
@ -214,7 +206,8 @@ subroutine CPFEM_init
call IO_read_realFile(777,'convergedTstar',modelName,size(crystallite_Tstar0_v)) call IO_read_realFile(777,'convergedTstar',modelName,size(crystallite_Tstar0_v))
read (777,rec=1) crystallite_Tstar0_v read (777,rec=1) crystallite_Tstar0_v
close (777) close (777)
#ifndef NEWSTATE
#ifdef TODO
call IO_read_realFile(777,'convergedStateConst',modelName) call IO_read_realFile(777,'convergedStateConst',modelName)
m = 0_pInt m = 0_pInt
do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems
@ -235,6 +228,8 @@ subroutine CPFEM_init
enddo; enddo enddo; enddo
close (777) close (777)
#endif #endif
call IO_read_realFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE)) call IO_read_realFile(777,'convergeddcsdE',modelName,size(CPFEM_dcsdE))
read (777,rec=1) CPFEM_dcsdE read (777,rec=1) CPFEM_dcsdE
close (777) close (777)
@ -303,19 +298,12 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
microstructure_elemhomo, & microstructure_elemhomo, &
#ifdef NEWSTATE plasticState, &
plasticState, & damageState, &
damageState, & thermalState, &
thermalState, & mappingConstitutive, &
mappingConstitutive, &
#endif
material_phase material_phase
#ifndef NEWSTATE
use constitutive, only: &
constitutive_state0, &
constitutive_state
#endif
use crystallite, only: & use crystallite, only: &
crystallite_partionedF,& crystallite_partionedF,&
crystallite_F0, & crystallite_F0, &
@ -402,14 +390,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
#ifndef NEWSTATE
forall ( i = 1:homogenization_maxNgrains, &
j = 1:mesh_maxNips, &
k = 1:mesh_NcpElems ) &
constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
#endif
#ifdef NEWSTATE
forall ( i = 1:size(plasticState)) plasticState(i)%state0= plasticState(i)%state ! copy state in this lenghty way because A component cannot be an array if the encompassing structure is an array forall ( i = 1:size(plasticState)) plasticState(i)%state0= plasticState(i)%state ! copy state in this lenghty way because A component cannot be an array if the encompassing structure is an array
forall ( i = 1:size(damageState)) damageState(i)%state0 = damageState(i)%state ! copy state in this lenghty way because A component cannot be an array if the encompassing structure is an array forall ( i = 1:size(damageState)) damageState(i)%state0 = damageState(i)%state ! copy state in this lenghty way because A component cannot be an array if the encompassing structure is an array
forall ( i = 1:size(thermalState)) thermalState(i)%state0= thermalState(i)%state ! copy state in this lenghty way because A component cannot be an array if the encompassing structure is an array forall ( i = 1:size(thermalState)) thermalState(i)%state0= thermalState(i)%state ! copy state in this lenghty way because A component cannot be an array if the encompassing structure is an array
@ -421,15 +401,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el
plasticState(mappingConstitutive(2,1,debug_i,debug_e))%state(:,mappingConstitutive(1,1,debug_i,debug_e)) plasticState(mappingConstitutive(2,1,debug_i,debug_e))%state(:,mappingConstitutive(1,1,debug_i,debug_e))
endif endif
endif endif
#else
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) then
write(6,'(a)') '<< CPFEM >> aging states'
if (debug_e <= mesh_NcpElems .and. debug_i <= mesh_maxNips) then
write(6,'(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)') '<< CPFEM >> aged state of elFE ip grain',&
debug_e, debug_i, 1, constitutive_state(1,debug_i,debug_e)%p
endif
endif
#endif
!$OMP PARALLEL DO !$OMP PARALLEL DO
do k = 1,mesh_NcpElems do k = 1,mesh_NcpElems
do j = 1,mesh_maxNips do j = 1,mesh_maxNips
@ -468,7 +440,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el
call IO_write_jobRealFile(777,'convergedTstar',size(crystallite_Tstar0_v)) call IO_write_jobRealFile(777,'convergedTstar',size(crystallite_Tstar0_v))
write (777,rec=1) crystallite_Tstar0_v write (777,rec=1) crystallite_Tstar0_v
close (777) close (777)
#ifndef NEWSTATE #ifdef TODO
call IO_write_jobRealFile(777,'convergedStateConst') call IO_write_jobRealFile(777,'convergedStateConst')
m = 0_pInt m = 0_pInt
do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems do i = 1,homogenization_maxNgrains; do j = 1,mesh_maxNips; do k = 1,mesh_NcpElems
@ -478,7 +450,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el
enddo enddo
enddo; enddo; enddo enddo; enddo; enddo
close (777) close (777)
#endif #endif
call IO_write_jobRealFile(777,'convergedStateHomog') call IO_write_jobRealFile(777,'convergedStateHomog')
m = 0_pInt m = 0_pInt
do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips do k = 1,mesh_NcpElems; do j = 1,mesh_maxNips
@ -648,8 +620,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature, dt, el
!*** copy to output if required (FEM solver) !*** copy to output if required (FEM solver)
if(present(cauchyStress)) cauchyStress = CPFEM_cs(1:6,ip,elCP) if(present(cauchyStress)) cauchyStress = CPFEM_cs (1:6, ip,elCP)
if(present(jacobian)) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) if(present(jacobian)) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP)
end subroutine CPFEM_general end subroutine CPFEM_general

35
code/Makefile
View File

@ -116,11 +116,6 @@ RUN_PATH :=$(RUN_PATH),-rpath,$(HDF5_ROOT)/lib
INCLUDE_DIRS +=-I$(HDF5_ROOT)/include -DHDF INCLUDE_DIRS +=-I$(HDF5_ROOT)/include -DHDF
endif endif
#new state
ifeq "$(STATE)" "NEW"
INCLUDE_DIRS +=-DNEWSTATE
endif
ifdef STANDARD_CHECK ifdef STANDARD_CHECK
STANDARD_CHECK_ifort =$(STANDARD_CHECK) STANDARD_CHECK_ifort =$(STANDARD_CHECK)
STANDARD_CHECK_gfortran =$(STANDARD_CHECK) STANDARD_CHECK_gfortran =$(STANDARD_CHECK)
@ -350,18 +345,22 @@ PRECISION_gfortran :=-fdefault-real-8 -fdefault-double-8 -DFLOAT=8 -DINT=4
COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) -DSpectral COMPILE =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(OPTI)_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) -DSpectral
COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) -DSpectral COMPILE_MAXOPTI =$(OPENMP_FLAG_$(F90)) $(COMPILE_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) $(INCLUDE_DIRS) $(PRECISION_$(F90)) -DSpectral
################################################################################################### ###################################################################################################
CONSTITUTIVE_FILES = constitutive_thermal.o \
constitutive_damage.o \
constitutive_nonlocal.o \
constitutive_titanmod.o \
constitutive_dislotwin.o \
constitutive_phenopowerlaw.o \
constitutive_j2.o \
constitutive_none.o
COMPILED_FILES = prec.o DAMASK_spectral_interface.o IO.o libs.o numerics.o debug.o math.o \ COMPILED_FILES = prec.o DAMASK_spectral_interface.o IO.o libs.o numerics.o debug.o math.o \
FEsolving.o mesh.o material.o lattice.o \ FEsolving.o mesh.o material.o lattice.o \
constitutive_dislotwin.o constitutive_j2.o constitutive_phenopowerlaw.o \ damage_none.o damage_gradient.o thermal_none.o thermal_conduction.o thermal_adiabatic.o \
constitutive_titanmod.o constitutive_nonlocal.o constitutive_none.o constitutive.o crystallite.o \ $(CONSTITUTIVE_FILES) constitutive.o crystallite.o \
homogenization_RGC.o homogenization_isostrain.o homogenization_none.o homogenization.o CPFEM.o \ homogenization_RGC.o homogenization_isostrain.o homogenization_none.o homogenization.o CPFEM.o \
DAMASK_spectral_utilities.o DAMASK_spectral_solverBasic.o DAMASK_spectral_utilities.o DAMASK_spectral_solverBasic.o \
ifeq "$(STATE)" "NEW"
COMPILED_FILES += constitutive_damage.o damage_none.o damage_gradient.o \
constitutive_thermal.o thermal_none.o thermal_conduction.o
CONSTITUTIVE_FILES := constitutive_thermal.o constitutive_damage.o
endif
ifdef PETSC_DIR ifdef PETSC_DIR
PETSC_FILES = DAMASK_spectral_solverAL.o \ PETSC_FILES = DAMASK_spectral_solverAL.o \
@ -410,15 +409,10 @@ homogenization_none.o: homogenization_none.f90 \
crystallite.o crystallite.o
crystallite.o: crystallite.f90 \ crystallite.o: crystallite.f90 \
constitutive.o $(CONSTITUTIVE_FILES) constitutive.o
constitutive.o: constitutive.f90 \ constitutive.o: constitutive.f90 \
constitutive_nonlocal.o \ $(CONSTITUTIVE_FILES)
constitutive_titanmod.o \
constitutive_dislotwin.o \
constitutive_phenopowerlaw.o \
constitutive_j2.o \
constitutive_none.o
constitutive_nonlocal.o: constitutive_nonlocal.f90 \ constitutive_nonlocal.o: constitutive_nonlocal.f90 \
lattice.o lattice.o
@ -450,6 +444,7 @@ damage_gradient.o: damage_gradient.f90 \
constitutive_thermal.o: constitutive_thermal.f90 \ constitutive_thermal.o: constitutive_thermal.f90 \
thermal_none.o \ thermal_none.o \
thermal_adiabatic.o \
thermal_conduction.o thermal_conduction.o
thermal_none.o: thermal_none.f90 \ thermal_none.o: thermal_none.f90 \

544
code/constitutive.f90
View File

@ -7,35 +7,12 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module constitutive module constitutive
use prec, only: & use prec, only: &
pInt, & pInt
pReal, &
p_vec
implicit none implicit none
private private
#ifndef NEWSTATE
type(p_vec), public, dimension(:,:,:), allocatable :: &
constitutive_state0, & !< pointer array to microstructure at start of BVP inc
constitutive_partionedState0, & !< pointer array to microstructure at start of homogenization inc
constitutive_subState0, & !< pointer array to microstructure at start of crystallite inc
constitutive_state, & !< pointer array to current microstructure (end of converged time step)
constitutive_state_backup, & !< pointer array to backed up microstructure (end of converged time step)
constitutive_dotState, & !< pointer array to evolution of current microstructure
constitutive_deltaState, & !< pointer array to incremental change of current microstructure
constitutive_previousDotState,& !< pointer array to previous evolution of current microstructure
constitutive_previousDotState2,& !< pointer array to 2nd previous evolution of current microstructure
constitutive_dotState_backup, & !< pointer array to backed up evolution of current microstructure
constitutive_RK4dotState, & !< pointer array to evolution of microstructure defined by classical Runge-Kutta method
constitutive_aTolState !< pointer array to absolute state tolerance
type(p_vec), public, dimension(:,:,:,:), allocatable :: &
constitutive_RKCK45dotState !< pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method
integer(pInt), public, dimension(:,:,:), allocatable :: & integer(pInt), public, dimension(:,:,:), allocatable :: &
constitutive_sizeDotState, & !< size of dotState array
constitutive_sizeState, & !< size of state array per grain
constitutive_sizePostResults !< size of postResults array per grain constitutive_sizePostResults !< size of postResults array per grain
integer(pInt), private :: &
constitutive_maxSizeState
#endif
integer(pInt), public, protected :: & integer(pInt), public, protected :: &
constitutive_maxSizePostResults, & constitutive_maxSizePostResults, &
constitutive_maxSizeDotState constitutive_maxSizeDotState
@ -68,6 +45,8 @@ subroutine constitutive_init
#endif #endif
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: &
pReal
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
@ -112,12 +91,9 @@ subroutine constitutive_init
PLASTICITY_PHENOPOWERLAW_label, & PLASTICITY_PHENOPOWERLAW_label, &
PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOTWIN_label, &
PLASTICITY_TITANMOD_label, & PLASTICITY_TITANMOD_label, &
#ifdef NEWSTATE
plasticState, & plasticState, &
#endif
#if defined(HDF) || defined(NEWSTATE)
mappingConstitutive, & mappingConstitutive, &
#endif
PLASTICITY_NONLOCAL_label PLASTICITY_NONLOCAL_label
use constitutive_none use constitutive_none
use constitutive_j2 use constitutive_j2
@ -214,32 +190,7 @@ subroutine constitutive_init
cMax = homogenization_maxNgrains cMax = homogenization_maxNgrains
iMax = mesh_maxNips iMax = mesh_maxNips
eMax = mesh_NcpElems eMax = mesh_NcpElems
#ifndef NEWSTATE allocate(constitutive_sizePostResults(cMax,iMax,eMax), source=0_pInt)
! lumped into new state
allocate(constitutive_state0(cMax,iMax,eMax))
allocate(constitutive_partionedState0(cMax,iMax,eMax))
allocate(constitutive_subState0(cMax,iMax,eMax))
allocate(constitutive_state(cMax,iMax,eMax))
allocate(constitutive_state_backup(cMax,iMax,eMax))
allocate(constitutive_dotState(cMax,iMax,eMax))
allocate(constitutive_deltaState(cMax,iMax,eMax))
allocate(constitutive_dotState_backup(cMax,iMax,eMax))
allocate(constitutive_aTolState(cMax,iMax,eMax))
! not needed anymore for new state
allocate(constitutive_sizeDotState(cMax,iMax,eMax), source=0_pInt)
allocate(constitutive_sizeState(cMax,iMax,eMax), source=0_pInt)
allocate(constitutive_sizePostResults(cMax,iMax,eMax), source=0_pInt)
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(cMax,iMax,eMax))
allocate(constitutive_previousDotState2(cMax,iMax,eMax))
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(cMax,iMax,eMax))
endif
if (any(numerics_integrator == 5_pInt)) then
allocate(constitutive_RKCK45dotState(6,cMax,iMax,eMax))
endif
ElemLoop:do e = 1_pInt,mesh_NcpElems ! loop over elements ElemLoop:do e = 1_pInt,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs IPloop:do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) ! loop over IPs
@ -251,207 +202,28 @@ subroutine constitutive_init
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
select case(phase_plasticity(material_phase(g,i,e))) select case(phase_plasticity(material_phase(g,i,e)))
case (PLASTICITY_NONE_ID) case (PLASTICITY_NONE_ID)
allocate(constitutive_state0(g,i,e)%p(constitutive_none_sizeState(instance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_none_sizeState(instance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_none_sizeState(instance)))
allocate(constitutive_state(g,i,e)%p(constitutive_none_sizeState(instance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_none_sizeState(instance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_none_sizeState(instance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_none_sizeDotState(instance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_none_sizeDotState(instance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_none_sizeDotState(instance)))
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_none_sizeDotState(instance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_none_sizeDotState(instance)))
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_none_sizeDotState(instance)))
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_none_sizeDotState(instance)))
enddo
endif
constitutive_state0(g,i,e)%p = 0.0_pReal
constitutive_aTolState(g,i,e)%p = 1.0_pReal
constitutive_sizeState(g,i,e) = 0_pInt
constitutive_sizeDotState(g,i,e) = 0_pInt
constitutive_sizePostResults(g,i,e) = 0_pInt constitutive_sizePostResults(g,i,e) = 0_pInt
case (PLASTICITY_J2_ID) case (PLASTICITY_J2_ID)
allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(instance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(instance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_j2_sizeState(instance)))
allocate(constitutive_state(g,i,e)%p(constitutive_j2_sizeState(instance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_j2_sizeState(instance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_j2_sizeState(instance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(instance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_j2_sizeDotState(instance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(instance)))
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(instance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(instance)))
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(instance)))
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(instance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(instance)
constitutive_aTolState(g,i,e)%p = constitutive_j2_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(instance)
constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(instance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(instance)
case (PLASTICITY_PHENOPOWERLAW_ID) case (PLASTICITY_PHENOPOWERLAW_ID)
allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_aTolState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_deltaState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(instance)),source=0.0_pReal)
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(instance)
constitutive_aTolState(g,i,e)%p = constitutive_phenopowerlaw_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(instance)
constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(instance) constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(instance)
case (PLASTICITY_DISLOTWIN_ID) case (PLASTICITY_DISLOTWIN_ID)
allocate(constitutive_state0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_subState0(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_state_backup(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_aTolState(g,i,e)%p(constitutive_dislotwin_sizeState(instance)),source=0.0_pReal)
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_deltaState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)),source=0.0_pReal)
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(instance,material_phase(g,i,e))
constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(instance)
constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(instance) constitutive_sizePostResults(g,i,e) = constitutive_dislotwin_sizePostResults(instance)
case (PLASTICITY_TITANMOD_ID) case (PLASTICITY_TITANMOD_ID)
allocate(constitutive_state0(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_state(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_titanmod_sizeState(instance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(instance,material_phase(g,i,e))
constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(instance)
constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(instance) constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(instance)
case (PLASTICITY_NONLOCAL_ID) case (PLASTICITY_NONLOCAL_ID)
nonlocalConstitutionPresent = .true. nonlocalConstitutionPresent = .true.
plasticState(mappingConstitutive(2,g,i,e))%nonlocal = .true.
if(myNgrains/=1_pInt) call IO_error(252_pInt, e,i,g) if(myNgrains/=1_pInt) call IO_error(252_pInt, e,i,g)
allocate(constitutive_state0(g,i,e)%p(constitutive_nonlocal_sizeState(instance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_nonlocal_sizeState(instance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_nonlocal_sizeState(instance)))
allocate(constitutive_state(g,i,e)%p(constitutive_nonlocal_sizeState(instance)))
allocate(constitutive_state_backup(g,i,e)%p(constitutive_nonlocal_sizeState(instance)))
allocate(constitutive_aTolState(g,i,e)%p(constitutive_nonlocal_sizeState(instance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(instance)))
allocate(constitutive_deltaState(g,i,e)%p(constitutive_nonlocal_sizeDotState(instance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(instance)))
if (any(numerics_integrator == 1_pInt)) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(instance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(instance)))
endif
if (any(numerics_integrator == 4_pInt)) then
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(instance)))
endif
if (any(numerics_integrator == 5_pInt)) then
do s = 1_pInt,6_pInt
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(instance)))
enddo
endif
constitutive_aTolState(g,i,e)%p = constitutive_nonlocal_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_nonlocal_sizeDotState(instance)
constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(instance) constitutive_sizePostResults(g,i,e) = constitutive_nonlocal_sizePostResults(instance)
end select end select
enddo GrainLoop enddo GrainLoop
enddo IPloop enddo IPloop
enddo ElemLoop enddo ElemLoop
#endif
#ifdef NEWSTATE
PhaseLoop:do phase = 1_pInt,material_Nphase ! loop over phases
instance = phase_plasticityInstance(phase)
select case(phase_plasticity(phase))
case (PLASTICITY_NONE_ID)
plasticState(phase)%sizePostResults = constitutive_none_sizePostResults(instance)
case (PLASTICITY_J2_ID)
plasticState(phase)%sizePostResults = constitutive_j2_sizePostResults(instance)
case (PLASTICITY_PHENOPOWERLAW_ID)
plasticState(phase)%sizePostResults = constitutive_none_sizePostResults(instance)
case (PLASTICITY_DISLOTWIN_ID)
plasticState(phase)%sizePostResults = constitutive_dislotwin_sizePostResults(instance)
case (PLASTICITY_TITANMOD_ID)
plasticState(phase)%sizePostResults = constitutive_titanmod_sizePostResults(instance)
case (PLASTICITY_NONLOCAL_ID)
nonlocalConstitutionPresent = .true.
plasticState(phase)%nonlocal = .true.
plasticState(phase)%sizePostResults = constitutive_nonlocal_sizePostResults(instance)
end select
enddo PhaseLoop
#endif
if (nonlocalConstitutionPresent) & if (nonlocalConstitutionPresent) &
#ifdef NEWSTATE
call constitutive_nonlocal_stateInit() call constitutive_nonlocal_stateInit()
#else
call constitutive_nonlocal_stateInit(constitutive_state0(1,1:iMax,1:eMax))
#endif
#ifdef NEWSTATE
do e = 1_pInt,mesh_NcpElems ! loop over elements do e = 1_pInt,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall(i = 1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))), g = 1_pInt:myNgrains) forall(i = 1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))), g = 1_pInt:myNgrains)
@ -461,16 +233,6 @@ subroutine constitutive_init
plasticState(mappingConstitutive(2,g,i,e))%State0(:,mappingConstitutive(1,g,i,e)) ! need to be defined for first call of constitutive_microstructure in crystallite_init plasticState(mappingConstitutive(2,g,i,e))%State0(:,mappingConstitutive(1,g,i,e)) ! need to be defined for first call of constitutive_microstructure in crystallite_init
endforall endforall
enddo enddo
#else
do e = 1_pInt,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall(i = 1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))), g = 1_pInt:myNgrains)
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p
constitutive_state(g,i,e)%p = constitutive_state0(g,i,e)%p ! need to be defined for first call of constitutive_microstructure in crystallite_init
endforall
enddo
#endif
#ifdef HDF #ifdef HDF
call HDF5_mappingConstitutive(mappingConstitutive) call HDF5_mappingConstitutive(mappingConstitutive)
@ -483,7 +245,7 @@ subroutine constitutive_init
enddo enddo
#endif #endif
#ifndef NEWSTATE #ifdef TODO
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! write out state size file ! write out state size file
call IO_write_jobIntFile(777,'sizeStateConst', size(constitutive_sizeState)) call IO_write_jobIntFile(777,'sizeStateConst', size(constitutive_sizeState))
@ -512,29 +274,30 @@ subroutine constitutive_init
write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', constitutive_maxSizePostResults write(6,'(a32,1x,7(i8,1x))') 'maxSizePostResults: ', constitutive_maxSizePostResults
endif endif
flush(6) flush(6)
#else #endif
constitutive_maxSizePostResults = 0_pInt constitutive_maxSizePostResults = 0_pInt
constitutive_maxSizeDotState = 0_pInt constitutive_maxSizeDotState = 0_pInt
do p = 1, size(plasticState) do p = 1, size(plasticState)
constitutive_maxSizeDotState = max(constitutive_maxSizeDotState, plasticState(p)%sizeDotState) constitutive_maxSizeDotState = max(constitutive_maxSizeDotState, plasticState(p)%sizeDotState)
constitutive_maxSizePostResults = max(constitutive_maxSizePostResults, plasticState(p)%sizePostResults) constitutive_maxSizePostResults = max(constitutive_maxSizePostResults, plasticState(p)%sizePostResults)
enddo enddo
#endif
end subroutine constitutive_init end subroutine constitutive_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_homogenizedC(ipc,ip,el) function constitutive_homogenizedC(ipc,ip,el)
use prec, only: &
pReal
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_phase, & material_phase, &
PLASTICITY_TITANMOD_ID, & PLASTICITY_TITANMOD_ID, &
#ifdef NEWSTATE
plasticState,& plasticState,&
mappingConstitutive, & mappingConstitutive, &
#endif
PLASTICITY_DISLOTWIN_ID PLASTICITY_DISLOTWIN_ID
use constitutive_titanmod, only: & use constitutive_titanmod, only: &
constitutive_titanmod_homogenizedC constitutive_titanmod_homogenizedC
@ -553,27 +316,11 @@ pure function constitutive_homogenizedC(ipc,ip,el)
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID) case (PLASTICITY_DISLOTWIN_ID)
#ifdef NEWSTATE constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(ipc,ip,el)
constitutive_homogenizedC = constitutive_dislotwin_homogenizedC &
(plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
constitutive_homogenizedC = constitutive_dislotwin_homogenizedC &
(constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_TITANMOD_ID) case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE constitutive_homogenizedC = constitutive_titanmod_homogenizedC (ipc,ip,el)
constitutive_homogenizedC = constitutive_titanmod_homogenizedC &
(plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state(ipc,ip,el), &
ipc,ip,el)
#endif
case default case default
constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phase(ipc,ip,el)) constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phase (ipc,ip,el))
end select end select
@ -584,14 +331,14 @@ end function constitutive_homogenizedC
!> @brief calls microstructure function of the different constitutive models !> @brief calls microstructure function of the different constitutive models
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el) subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
use prec, only: &
pReal
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_phase, & material_phase, &
PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOTWIN_ID, &
#ifdef NEWSTATE
plasticState, & plasticState, &
mappingConstitutive, & mappingConstitutive, &
#endif
PLASTICITY_TITANMOD_ID, & PLASTICITY_TITANMOD_ID, &
PLASTICITY_NONLOCAL_ID PLASTICITY_NONLOCAL_ID
use constitutive_titanmod, only: & use constitutive_titanmod, only: &
@ -615,30 +362,12 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el)
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_DISLOTWIN_ID) case (PLASTICITY_DISLOTWIN_ID)
#ifdef NEWSTATE call constitutive_dislotwin_microstructure(temperature,ipc,ip,el)
call constitutive_dislotwin_microstructure(temperature, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
call constitutive_dislotwin_microstructure(temperature,constitutive_state(ipc,ip,el), &
ipc,ip,el)
#endif
case (PLASTICITY_TITANMOD_ID) case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE call constitutive_titanmod_microstructure (temperature,ipc,ip,el)
call constitutive_titanmod_microstructure(temperature, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
call constitutive_titanmod_microstructure(temperature,constitutive_state(ipc,ip,el), &
ipc,ip,el)
#endif
case (PLASTICITY_NONLOCAL_ID) case (PLASTICITY_NONLOCAL_ID)
#ifdef NEWSTATE call constitutive_nonlocal_microstructure (Fe,Fp, ip,el)
call constitutive_nonlocal_microstructure(Fe,Fp,ipc,ip,el)
#else
call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el)
#endif
end select end select
end subroutine constitutive_microstructure end subroutine constitutive_microstructure
@ -648,15 +377,15 @@ end subroutine constitutive_microstructure
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ipc, ip, el) subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ipc, ip, el)
use prec, only: &
pReal
use math, only: & use math, only: &
math_identity2nd math_identity2nd
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
material_phase, & material_phase, &
#ifdef NEWSTATE
plasticState,& plasticState,&
mappingConstitutive, & mappingConstitutive, &
#endif
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
PLASTICITY_J2_ID, & PLASTICITY_J2_ID, &
PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_PHENOPOWERLAW_ID, &
@ -693,50 +422,17 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip
case (PLASTICITY_NONE_ID) case (PLASTICITY_NONE_ID)
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar = math_identity2nd(9) dLp_dTstar = math_identity2nd(9)
case (PLASTICITY_J2_ID) case (PLASTICITY_J2_ID)
#ifdef NEWSTATE call constitutive_j2_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,ipc,ip,el)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)),ipc,ip,el)
#else
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_PHENOPOWERLAW_ID) case (PLASTICITY_PHENOPOWERLAW_ID)
#ifdef NEWSTATE call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,ipc,ip,el)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)),ipc,ip,el)
#else
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_DISLOTWIN_ID)
#ifdef NEWSTATE
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
temperature,constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,temperature, &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(:,mappingConstitutive(1,ipc,ip,el)), &
ipc,ip,el)
#else
call constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, &
temperature,constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_NONLOCAL_ID) case (PLASTICITY_NONLOCAL_ID)
#ifdef NEWSTATE call constitutive_nonlocal_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature, ip,el)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, & case (PLASTICITY_DISLOTWIN_ID)
temperature, ipc,ip,el) call constitutive_dislotwin_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
#else case (PLASTICITY_TITANMOD_ID)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, & call constitutive_titanmod_LpAndItsTangent (Lp,dLp_dTstar,Tstar_v,temperature,ipc,ip,el)
temperature, constitutive_state(ipc,ip,el), ipc,ip,el)
#endif
end select end select
end subroutine constitutive_LpAndItsTangent end subroutine constitutive_LpAndItsTangent
@ -748,8 +444,10 @@ end subroutine constitutive_LpAndItsTangent
!> the elastic deformation gradient depending on the selected elastic law (so far no case switch !> the elastic deformation gradient depending on the selected elastic law (so far no case switch
!! because only hooke is implemented !! because only hooke is implemented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) subroutine constitutive_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
use prec, only: &
pReal
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
@ -772,28 +470,28 @@ end subroutine constitutive_TandItsTangent
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic deformation gradient using hookes law !> the elastic deformation gradient using hookes law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el) subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, Fe, ipc, ip, el)
use math, only : & use prec, only: &
math_mul3x3, & pReal
math_mul33x33, & use math, only : &
math_mul3333xx33, & math_mul3x3, &
math_Mandel66to3333, & math_mul33x33, &
math_transpose33, & math_mul3333xx33, &
MATH_I3 math_Mandel66to3333, &
#ifdef NEWSTATE math_transpose33, &
use material, only: & MATH_I3
mappingConstitutive, & use material, only: &
damageState, & mappingConstitutive, &
phase_damage, & damageState, &
DAMAGE_gradient_ID, & phase_damage, &
thermalState, & DAMAGE_gradient_ID, &
phase_thermal, & thermalState, &
THERMAL_conduction_ID, & phase_thermal, &
THERMAL_adiabatic_ID THERMAL_conduction_ID, &
THERMAL_adiabatic_ID
use lattice, only: & use lattice, only: &
lattice_referenceTemperature, & lattice_referenceTemperature, &
lattice_thermalExpansion33 lattice_thermalExpansion33
#endif
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -822,7 +520,6 @@ use material, only: &
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) & forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt, k=1_pInt:3_pInt, l=1_pInt:3_pInt) &
dT_dFe(i,j,k,l) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3)) ! dT*_ij/dFe_kl dT_dFe(i,j,k,l) = math_mul3x3(C(i,j,l,1:3),Fe(k,1:3)) ! dT*_ij/dFe_kl
#ifdef NEWSTATE
phase = mappingConstitutive(2,ipc,ip,el) phase = mappingConstitutive(2,ipc,ip,el)
constituent = mappingConstitutive(1,ipc,ip,el) constituent = mappingConstitutive(1,ipc,ip,el)
select case (phase_damage(phase)) select case (phase_damage(phase))
@ -842,7 +539,6 @@ use material, only: &
lattice_referenceTemperature(phase)) & lattice_referenceTemperature(phase)) &
* lattice_thermalExpansion33(1:3,1:3,phase)) * lattice_thermalExpansion33(1:3,1:3,phase))
end select end select
#endif
end subroutine constitutive_hooke_TandItsTangent end subroutine constitutive_hooke_TandItsTangent
@ -853,6 +549,7 @@ end subroutine constitutive_hooke_TandItsTangent
subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, subdt, subfracArray,& subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, subdt, subfracArray,&
ipc, ip, el) ipc, ip, el)
use prec, only: & use prec, only: &
pReal, &
pLongInt pLongInt
use debug, only: & use debug, only: &
debug_cumDotStateCalls, & debug_cumDotStateCalls, &
@ -865,10 +562,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
phase_plasticity, & phase_plasticity, &
#ifdef NEWSTATE
plasticState, & plasticState, &
mappingConstitutive, & mappingConstitutive, &
#endif
material_phase, & material_phase, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
@ -912,59 +607,17 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature,
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID)
case (PLASTICITY_J2_ID) case (PLASTICITY_J2_ID)
#ifdef NEWSTATE call constitutive_j2_dotState (Tstar_v,ipc,ip,el)
plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) &
= constitutive_j2_dotState(Tstar_v,plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)), ipc,ip,el)
#else
constitutive_dotState(ipc,ip,el)%p = constitutive_j2_dotState(Tstar_v,&
constitutive_state(ipc,ip,el), ipc,ip,el)
#endif
case (PLASTICITY_PHENOPOWERLAW_ID) case (PLASTICITY_PHENOPOWERLAW_ID)
#ifdef NEWSTATE call constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) &
= constitutive_phenopowerlaw_dotState(Tstar_v,plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)), ipc,ip,el)
#else
constitutive_dotState(ipc,ip,el)%p = constitutive_phenopowerlaw_dotState(Tstar_v,&
constitutive_state(ipc,ip,el), ipc,ip,el)
#endif
case (PLASTICITY_DISLOTWIN_ID) case (PLASTICITY_DISLOTWIN_ID)
#ifdef NEWSTATE call constitutive_dislotwin_dotState (Tstar_v,Temperature,ipc,ip,el)
plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) &
= constitutive_dislotwin_dotState(Tstar_v,Temperature,&
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)), ipc,ip,el)
#else
constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,&
constitutive_state(ipc,ip,el), ipc,ip,el)
#endif
case (PLASTICITY_TITANMOD_ID) case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE call constitutive_titanmod_dotState (Tstar_v,Temperature,ipc,ip,el)
plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) &
= constitutive_titanmod_dotState(Tstar_v,Temperature,&
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)), ipc,ip,el)
#else
constitutive_dotState(ipc,ip,el)%p = constitutive_titanmod_dotState(Tstar_v,Temperature,&
constitutive_state(ipc,ip,el), ipc,ip,el)
#endif
case (PLASTICITY_NONLOCAL_ID) case (PLASTICITY_NONLOCAL_ID)
#ifdef NEWSTATE call constitutive_nonlocal_dotState (Tstar_v,FeArray,FpArray,Temperature, subdt, &
!* subfracArray,ip,el)
plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) = &
constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, &
Temperature, subdt, &
subfracArray, ipc, ip, el)
#else
constitutive_dotState(ipc,ip,el)%p = constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, &
Temperature, constitutive_state, constitutive_state0, subdt, &
subfracArray, ipc, ip, el)
#endif
end select end select
if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then
@ -984,6 +637,7 @@ end subroutine constitutive_collectDotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
logical function constitutive_collectDeltaState(Tstar_v, ipc, ip, el) logical function constitutive_collectDeltaState(Tstar_v, ipc, ip, el)
use prec, only: & use prec, only: &
pReal, &
pLongInt pLongInt
use debug, only: & use debug, only: &
debug_cumDeltaStateCalls, & debug_cumDeltaStateCalls, &
@ -1021,17 +675,10 @@ logical function constitutive_collectDeltaState(Tstar_v, ipc, ip, el)
case (PLASTICITY_NONLOCAL_ID) case (PLASTICITY_NONLOCAL_ID)
constitutive_collectDeltaState = .true. constitutive_collectDeltaState = .true.
#ifdef NEWSTATE
call constitutive_nonlocal_deltaState(Tstar_v,ip,el) call constitutive_nonlocal_deltaState(Tstar_v,ip,el)
#else
call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el)%p,&
constitutive_state(ipc,ip,el), Tstar_v,ipc,ip,el)
#endif
case default case default
constitutive_collectDeltaState = .false. constitutive_collectDeltaState = .false.
#ifndef NEWSTATE
constitutive_deltaState(ipc,ip,el)%p = 0.0_pReal !ToDo: needed or will it remain zero anyway?
#endif
end select end select
if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then
@ -1051,14 +698,14 @@ end function constitutive_collectDeltaState
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el) function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
use prec, only: &
pReal
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
#ifdef NEWSTATE
plasticState, & plasticState, &
mappingConstitutive, & mappingConstitutive, &
#endif
phase_plasticity, & phase_plasticity, &
material_phase, & material_phase, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
@ -1086,13 +733,8 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
ipc, & !< grain number ipc, & !< grain number
ip, & !< integration point number ip, & !< integration point number
el !< element number el !< element number
#ifndef NEWSTATE
real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: & real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: &
constitutive_postResults constitutive_postResults
#else
real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults) :: &
constitutive_postResults
#endif
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
temperature temperature
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
@ -1103,52 +745,16 @@ function constitutive_postResults(Tstar_v, FeArray, temperature, ipc, ip, el)
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
select case (phase_plasticity(material_phase(ipc,ip,el))) select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID)
case (PLASTICITY_TITANMOD_ID) case (PLASTICITY_TITANMOD_ID)
#ifdef NEWSTATE constitutive_postResults = constitutive_titanmod_postResults (ipc,ip,el)
constitutive_postResults = constitutive_titanmod_postResults(&
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)),ipc,ip,el)
#else
constitutive_postResults = constitutive_titanmod_postResults(&
constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_J2_ID) case (PLASTICITY_J2_ID)
#ifdef NEWSTATE constitutive_postResults= constitutive_j2_postResults (Tstar_v,ipc,ip,el)
constitutive_postResults= constitutive_j2_postResults(Tstar_v, &
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)),ipc,ip,el)
#else
constitutive_postResults = constitutive_j2_postResults(Tstar_v,&
constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_PHENOPOWERLAW_ID) case (PLASTICITY_PHENOPOWERLAW_ID)
#ifdef NEWSTATE constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,&
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)),ipc,ip,el)
#else
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,&
constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_DISLOTWIN_ID) case (PLASTICITY_DISLOTWIN_ID)
#ifdef NEWSTATE constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,&
plasticState(mappingConstitutive(2,ipc,ip,el))% &
state(:,mappingConstitutive(1,ipc,ip,el)),ipc,ip,el)
#else
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,&
constitutive_state(ipc,ip,el),ipc,ip,el)
#endif
case (PLASTICITY_NONLOCAL_ID) case (PLASTICITY_NONLOCAL_ID)
#ifdef NEWSTATE constitutive_postResults = constitutive_nonlocal_postResults (Tstar_v,FeArray, ip,el)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, FeArray, &
mappingConstitutive, ipc, ip, el)
#else
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, FeArray, &
constitutive_state, constitutive_dotstate(ipc,ip,el), ipc, ip, el)
#endif
end select end select
end function constitutive_postResults end function constitutive_postResults

23
code/constitutive_damage.f90
View File

@ -19,7 +19,6 @@ module constitutive_damage
constitutive_damage_init, & constitutive_damage_init, &
constitutive_damage_microstructure, & constitutive_damage_microstructure, &
constitutive_damage_collectDotState, & constitutive_damage_collectDotState, &
constitutive_damage_collectDeltaState, &
constitutive_damage_postResults constitutive_damage_postResults
contains contains
@ -206,28 +205,6 @@ subroutine constitutive_damage_collectDotState(Tstar_v, Lp, ipc, ip, el)
end subroutine constitutive_damage_collectDotState end subroutine constitutive_damage_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state (so far, only nonlocal)
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
logical function constitutive_damage_collectDeltaState(ipc, ip, el)
use material, only: &
material_phase, &
phase_damage
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
select case (phase_damage(material_phase(ipc,ip,el)))
end select
constitutive_damage_collectDeltaState = .true.
end function constitutive_damage_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results

679
code/constitutive_dislotwin.f90
View File

File diff suppressed because it is too large Load Diff

225
code/constitutive_j2.f90
View File

@ -13,6 +13,7 @@ module constitutive_j2
use hdf5, only: & use hdf5, only: &
HID_T HID_T
#endif #endif
use prec, only: & use prec, only: &
pReal,& pReal,&
pInt pInt
@ -20,10 +21,6 @@ module constitutive_j2
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
#ifndef NEWSTATE
constitutive_j2_sizeDotState, & !< number of dotStates
constitutive_j2_sizeState, & !< total number of microstructural variables
#endif
constitutive_j2_sizePostResults !< cumulative size of post results constitutive_j2_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -58,25 +55,24 @@ module constitutive_j2
flowstress_ID, & flowstress_ID, &
strainrate_ID strainrate_ID
end enum end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
constitutive_j2_outputID !< ID of each post result output constitutive_j2_outputID !< ID of each post result output
#ifdef HDF #ifdef HDF
type constitutive_j2_tOutput type constitutive_j2_tOutput
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
flowstress, & flowstress, &
strainrate 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)) 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 constitutive_j2_tOutput end type constitutive_j2_tOutput
type(constitutive_j2_tOutput), allocatable, dimension(:) :: constitutive_j2_Output2 type(constitutive_j2_tOutput), allocatable, dimension(:) :: constitutive_j2_Output2
integer(HID_T), allocatable, dimension(:) :: outID integer(HID_T), allocatable, dimension(:) :: outID
#endif #endif
public :: & public :: &
constitutive_j2_init, & constitutive_j2_init, &
#ifndef NEWSTATE
constitutive_j2_stateInit, &
constitutive_j2_aTolState, &
#endif
constitutive_j2_LpAndItsTangent, & constitutive_j2_LpAndItsTangent, &
constitutive_j2_dotState, & constitutive_j2_dotState, &
constitutive_j2_postResults constitutive_j2_postResults
@ -126,9 +122,7 @@ subroutine constitutive_j2_init(fileUnit)
PLASTICITY_J2_label, & PLASTICITY_J2_label, &
PLASTICITY_J2_ID, & PLASTICITY_J2_ID, &
material_phase, & material_phase, &
#ifdef NEWSTATE
plasticState, & plasticState, &
#endif
MATERIAL_partPhase MATERIAL_partPhase
use lattice use lattice
@ -139,19 +133,25 @@ subroutine constitutive_j2_init(fileUnit)
integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt) :: phase, maxNinstance, instance,o, mySize, myConstituents integer(pInt) :: &
o, &
phase, &
maxNinstance, &
instance, &
mySize, &
sizeDotState, &
sizeState
character(len=65536) :: & character(len=65536) :: &
tag = '', & tag = '', &
line = '' line = ''
integer(pInt) :: NofMyPhase integer(pInt) :: NofMyPhase
#ifdef HDF #ifdef HDF
character(len=5) :: & character(len=5) :: &
str1 str1
integer(HID_T) :: ID,ID2,ID4 integer(HID_T) :: ID,ID2,ID4
#endif #endif
#ifdef NEWSTATE
integer(pInt) :: sizeDotState,sizeState
#endif
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_J2_label//' init -+>>>' write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_J2_label//' init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -160,17 +160,14 @@ subroutine constitutive_j2_init(fileUnit)
maxNinstance = int(count(phase_plasticity == PLASTICITY_J2_ID),pInt) maxNinstance = int(count(phase_plasticity == PLASTICITY_J2_ID),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifdef HDF #ifdef HDF
allocate(constitutive_j2_Output2(maxNinstance)) allocate(constitutive_j2_Output2(maxNinstance))
allocate(outID(maxNinstance)) allocate(outID(maxNinstance))
#endif #endif
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifndef NEWSTATE
allocate(constitutive_j2_sizeDotState(maxNinstance), source=1_pInt)
allocate(constitutive_j2_sizeState(maxNinstance), source=1_pInt)
#endif
allocate(constitutive_j2_sizePostResults(maxNinstance), source=0_pInt) allocate(constitutive_j2_sizePostResults(maxNinstance), source=0_pInt)
allocate(constitutive_j2_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt) allocate(constitutive_j2_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
allocate(constitutive_j2_output(maxval(phase_Noutput), maxNinstance)) allocate(constitutive_j2_output(maxval(phase_Noutput), maxNinstance))
@ -205,18 +202,16 @@ subroutine constitutive_j2_init(fileUnit)
exit exit
endif endif
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
myConstituents = 0_pInt
phase = phase + 1_pInt ! advance section counter phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_J2_ID) then if (phase_plasticity(phase) == PLASTICITY_J2_ID) then
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
myConstituents = count(material_phase==phase)
#ifdef HDF #ifdef HDF
outID(instance)=HDF5_addGroup(str1,tempResults) outID(instance)=HDF5_addGroup(str1,tempResults)
#endif #endif
endif endif
cycle ! skip to next line cycle ! skip to next line
endif endif
if (myConstituents > 0_pInt ) then 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 instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
positions = IO_stringPos(line,MAXNCHUNKS) positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
@ -290,13 +285,15 @@ subroutine constitutive_j2_init(fileUnit)
case default case default
end select end select
endif endif; endif
enddo parsingFile enddo parsingFile
initializeInstances: do phase = 1_pInt, size(phase_plasticity) initializeInstances: do phase = 1_pInt, size(phase_plasticity)
NofMyPhase=count(material_phase==phase) myPhase: if (phase_plasticity(phase) == PLASTICITY_j2_ID) then
if (phase_plasticity(phase) == PLASTICITY_j2_ID .and. NofMyPhase/=0) then NofMyPhase=count(material_phase==phase)
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,constitutive_j2_Noutput(instance) outputsLoop: do o = 1_pInt,constitutive_j2_Noutput(instance)
select case(constitutive_j2_outputID(o,instance)) select case(constitutive_j2_outputID(o,instance))
case(flowstress_ID,strainrate_ID) case(flowstress_ID,strainrate_ID)
@ -310,75 +307,41 @@ subroutine constitutive_j2_init(fileUnit)
constitutive_j2_sizePostResults(instance) + mySize constitutive_j2_sizePostResults(instance) + mySize
endif endif
enddo outputsLoop enddo outputsLoop
#ifdef NEWSTATE
sizeState = 1 !--------------------------------------------------------------------------------------------------
! allocate state arrays
sizeState = 1_pInt
plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeState = sizeState
sizeDotState = sizeState sizeDotState = sizeState
plasticState(phase)%sizeDotState = sizeDotState plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizePostResults = constitutive_j2_sizePostResults(instance) plasticState(phase)%sizePostResults = constitutive_j2_sizePostResults(instance)
allocate(plasticState(phase)%state0 (sizeState,NofMyPhase),source=constitutive_j2_tau0(instance)) allocate(plasticState(phase)%aTolState (sizeState),source=constitutive_j2_aTolResistance(instance))
allocate(plasticState(phase)%partionedState0(sizeState,NofMyPhase),source=constitutive_j2_tau0(instance)) allocate(plasticState(phase)%state0 (sizeState,NofMyPhase),source=constitutive_j2_tau0(instance))
allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%state (sizeState,NofMyPhase),source=constitutive_j2_tau0(instance)) allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%state (sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%aTolState (NofMyPhase),source=constitutive_j2_aTolResistance(instance)) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase),source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase),source=0.0_pReal)
endif endif
if (any(numerics_integrator == 4_pInt)) & if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase),source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) & if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
#endif
endif endif myPhase
enddo initializeInstances enddo initializeInstances
end subroutine constitutive_j2_init end subroutine constitutive_j2_init
#ifndef NEWSTATE
!--------------------------------------------------------------------------------------------------
!> @brief sets the initial microstructural state for a given instance of this plasticity
!> @details initial microstructural state is set to the value specified by tau0
! not needed for new state
!--------------------------------------------------------------------------------------------------
pure function constitutive_j2_stateInit(instance)
implicit none
real(pReal), dimension(1) :: constitutive_j2_stateInit
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
constitutive_j2_stateInit = constitutive_j2_tau0(instance)
end function constitutive_j2_stateInit
!--------------------------------------------------------------------------------------------------
!> @brief sets the relevant state values for a given instance of this plasticity
! not needed for new state
!--------------------------------------------------------------------------------------------------
pure function constitutive_j2_aTolState(instance)
implicit none
real(pReal), dimension(1) :: constitutive_j2_aTolState
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
constitutive_j2_aTolState = constitutive_j2_aTolResistance(instance)
end function constitutive_j2_aTolState
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ipc,ip,el) subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use prec, only: &
p_vec
use math, only: & use math, only: &
math_mul6x6, & math_mul6x6, &
math_Mandel6to33, & math_Mandel6to33, &
@ -389,6 +352,8 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ip
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
mappingConstitutive, &
plasticState, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
@ -405,13 +370,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ip
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
#ifdef NEWSTATE
real(pReal), dimension(1), intent(in) :: &
state
#else
type(p_vec), intent(in) :: &
state !< microstructure state
#endif
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
@ -434,15 +392,11 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ip
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal dLp_dTstar99 = 0.0_pReal
else else
#ifdef NEWSTATE
gamma_dot = constitutive_j2_gdot0(instance) & gamma_dot = constitutive_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / constitutive_j2_fTaylor(instance) / state(1)) & * (sqrt(1.5_pReal) * norm_Tstar_dev / (constitutive_j2_fTaylor(instance) * &
plasticState(mappingConstitutive(2,ipc,ip,el))%state(1,mappingConstitutive(1,ipc,ip,el)))) &
**constitutive_j2_n(instance) **constitutive_j2_n(instance)
#else
gamma_dot = constitutive_j2_gdot0(instance) &
* (sqrt(1.5_pReal) * norm_Tstar_dev / constitutive_j2_fTaylor(instance) / state%p(1)) &
**constitutive_j2_n(instance)
#endif
Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(instance) Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(instance)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -463,23 +417,21 @@ end subroutine constitutive_j2_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_dotState(Tstar_v,state,ipc,ip,el) subroutine constitutive_j2_dotState(Tstar_v,ipc,ip,el)
use prec, only: &
p_vec
use math, only: & use math, only: &
math_mul6x6 math_mul6x6
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
mappingConstitutive, &
plasticState, &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(1) :: & real(pReal) :: &
constitutive_j2_dotState
real(pReal) :: &
tempState tempState
real(pReal), dimension(6), intent(in):: & real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
@ -487,13 +439,6 @@ real(pReal) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
#ifdef NEWSTATE
real(pReal), dimension(1), intent(in) :: &
state
#else
type(p_vec), intent(in) :: &
state !< microstructure state
#endif
real(pReal), dimension(6) :: & real(pReal), dimension(6) :: &
Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal) :: & real(pReal) :: &
@ -502,14 +447,14 @@ real(pReal) :: &
saturation, & !< saturation resistance saturation, & !< saturation resistance
norm_Tstar_dev !< euclidean norm of Tstar_dev norm_Tstar_dev !< euclidean norm of Tstar_dev
integer(pInt) :: & integer(pInt) :: &
instance instance, & !< instance of my instance (unique number of my constitutive model)
#ifdef NEWSTATE of, & !< shortcut notation for offset position in state array
tempState = state(1) ph !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
#else
tempState = state%p(1)
#endif
of = mappingConstitutive(1,ipc,ip,el)
ph = mappingConstitutive(2,ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! norm of deviatoric part of 2nd Piola-Kirchhoff stress ! 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(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
@ -518,9 +463,9 @@ real(pReal) :: &
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! strain rate ! strain rate
gamma_dot = constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & gamma_dot = constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!----------------------------------------------------------------------------------- / &!-----------------------------------------------------------------------------------
(constitutive_j2_fTaylor(instance) * tempState) ) ** constitutive_j2_n(instance) (constitutive_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**constitutive_j2_n(instance)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! hardening coefficient ! hardening coefficient
@ -543,21 +488,20 @@ real(pReal) :: &
endif endif
hardening = ( constitutive_j2_h0(instance) + constitutive_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) & hardening = ( constitutive_j2_h0(instance) + constitutive_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) &
* abs( 1.0_pReal - tempState/saturation )**constitutive_j2_a(instance) & * abs( 1.0_pReal - tempState/saturation )**constitutive_j2_a(instance) &
* sign(1.0_pReal, 1.0_pReal - tempState/saturation) * sign(1.0_pReal, 1.0_pReal - plasticState(ph)%state(1,of)/saturation)
else else
hardening = 0.0_pReal hardening = 0.0_pReal
endif endif
plasticState(ph)%dotState(1,of) = hardening * gamma_dot !!!!!!!!!!!!!check if dostate
end subroutine constitutive_j2_dotState
constitutive_j2_dotState = hardening * gamma_dot
end function constitutive_j2_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el) function constitutive_j2_postResults(Tstar_v,ipc,ip,el)
use prec, only: &
p_vec
use math, only: & use math, only: &
math_mul6x6 math_mul6x6
use mesh, only: & use mesh, only: &
@ -566,6 +510,8 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el)
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase, & material_phase, &
plasticState, &
mappingConstitutive, &
phase_plasticityInstance, & phase_plasticityInstance, &
phase_Noutput phase_Noutput
@ -576,15 +522,6 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal) :: &
tempState
#ifdef NEWSTATE
real(pReal), dimension(1), intent(in) :: &
state
#else
type(p_vec), intent(in) :: &
state !< microstructure state
#endif
real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_j2_postResults constitutive_j2_postResults
@ -593,16 +530,14 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el)
real(pReal) :: & real(pReal) :: &
norm_Tstar_dev ! euclidean norm of Tstar_dev norm_Tstar_dev ! euclidean norm of Tstar_dev
integer(pInt) :: & integer(pInt) :: &
instance, & instance, & !< instance of my instance (unique number of my constitutive model)
o, & of, & !< shortcut notation for offset position in state array
c ph, & !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
c, &
o
#ifdef NEWSTATE of = mappingConstitutive(1,ipc,ip,el)
tempState = state(1) ph = mappingConstitutive(2,ipc,ip,el)
#else
tempState = state%p(1)
#endif
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -617,13 +552,13 @@ pure function constitutive_j2_postResults(Tstar_v,state,ipc,ip,el)
outputsLoop: do o = 1_pInt,constitutive_j2_Noutput(instance) outputsLoop: do o = 1_pInt,constitutive_j2_Noutput(instance)
select case(constitutive_j2_outputID(o,instance)) select case(constitutive_j2_outputID(o,instance))
case (flowstress_ID) case (flowstress_ID)
constitutive_j2_postResults(c+1_pInt) = tempState constitutive_j2_postResults(c+1_pInt) = plasticState(ph)%state(1,of)
c = c + 1_pInt c = c + 1_pInt
case (strainrate_ID) case (strainrate_ID)
constitutive_j2_postResults(c+1_pInt) = & constitutive_j2_postResults(c+1_pInt) = &
constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & constitutive_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev &
/ &!---------------------------------------------------------------------------------- / &!----------------------------------------------------------------------------------
(constitutive_j2_fTaylor(instance) * tempState) ) ** constitutive_j2_n(instance) (constitutive_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** constitutive_j2_n(instance)
c = c + 1_pInt c = c + 1_pInt
end select end select
enddo outputsLoop enddo outputsLoop

47
code/constitutive_none.f90
View File

@ -12,10 +12,6 @@ module constitutive_none
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
#ifndef NEWSTATE
constitutive_none_sizeDotState, &
constitutive_none_sizeState, &
#endif
constitutive_none_sizePostResults constitutive_none_sizePostResults
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -45,11 +41,8 @@ subroutine constitutive_none_init(fileUnit)
phase_plasticity, & phase_plasticity, &
phase_Noutput, & phase_Noutput, &
PLASTICITY_NONE_label, & PLASTICITY_NONE_label, &
#ifdef NEWSTATE
material_phase, & material_phase, &
plasticState, & plasticState, &
phase_plasticityInstance, &
#endif
PLASTICITY_none_ID, & PLASTICITY_none_ID, &
MATERIAL_partPhase MATERIAL_partPhase
@ -57,7 +50,6 @@ subroutine constitutive_none_init(fileUnit)
integer(pInt), intent(in) :: fileUnit integer(pInt), intent(in) :: fileUnit
integer(pInt) :: & integer(pInt) :: &
instance, &
maxNinstance, & maxNinstance, &
phase, & phase, &
NofMyPhase, & NofMyPhase, &
@ -75,39 +67,36 @@ subroutine constitutive_none_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(constitutive_none_sizePostResults(maxNinstance), source=0_pInt)
#ifdef NEWSTATE
initializeInstances: do phase = 1_pInt, size(phase_plasticity) initializeInstances: do phase = 1_pInt, size(phase_plasticity)
if (phase_plasticity(phase) == PLASTICITY_none_ID) then
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
if (phase_plasticity(phase) == PLASTICITY_none_ID .and. NofMyPhase/=0) then
instance = phase_plasticityInstance(phase)
sizeState = 0_pInt sizeState = 0_pInt
plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeState = sizeState
sizeDotState = sizeState sizeDotState = sizeState
plasticState(phase)%sizeDotState = sizeDotState plasticState(phase)%sizeDotState = sizeDotState
plasticState(phase)%sizePostResults = constitutive_none_sizePostResults(instance) plasticState(phase)%sizePostResults = 0_pInt
allocate(plasticState(phase)%state0 (sizeState,NofMyPhase)) allocate(plasticState(phase)%aTolState (sizeState))
allocate(plasticState(phase)%partionedState0(sizeState,NofMyPhase)) allocate(plasticState(phase)%state0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase)) allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%state (sizeState,NofMyPhase)) allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase))
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase)) allocate(plasticState(phase)%state (sizeState,NofMyPhase))
allocate(plasticState(phase)%aTolState (NofMyPhase)) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase))
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase))
allocate(plasticState(phase)%dotState_backup(sizeDotState,NofMyPhase)) allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase))
allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase))
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase)) allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase))
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase)) allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase))
endif endif
if (any(numerics_integrator == 4_pInt)) & if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase)) allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase))
if (any(numerics_integrator == 5_pInt)) & if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase)) allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase))
endif endif
enddo initializeInstances enddo initializeInstances
#else
allocate(constitutive_none_sizeDotState(maxNinstance), source=1_pInt) allocate(constitutive_none_sizePostResults(maxNinstance), source=0_pInt)
allocate(constitutive_none_sizeState(maxNinstance), source=1_pInt)
#endif
end subroutine constitutive_none_init end subroutine constitutive_none_init

1519
code/constitutive_nonlocal.f90
View File

File diff suppressed because it is too large Load Diff

539
code/constitutive_phenopowerlaw.f90
View File

@ -14,10 +14,6 @@ module constitutive_phenopowerlaw
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
#ifndef NEWSTATE
constitutive_phenopowerlaw_sizeDotState, &
constitutive_phenopowerlaw_sizeState, &
#endif
constitutive_phenopowerlaw_sizePostResults !< cumulative size of post results constitutive_phenopowerlaw_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -40,7 +36,6 @@ module constitutive_phenopowerlaw
constitutive_phenopowerlaw_gdot0_twin, & !< reference shear strain rate for twin (input parameter) constitutive_phenopowerlaw_gdot0_twin, & !< reference shear strain rate for twin (input parameter)
constitutive_phenopowerlaw_n_slip, & !< stress exponent for slip (input parameter) constitutive_phenopowerlaw_n_slip, & !< stress exponent for slip (input parameter)
constitutive_phenopowerlaw_n_twin, & !< stress exponent for twin (input parameter) constitutive_phenopowerlaw_n_twin, & !< stress exponent for twin (input parameter)
constitutive_phenopowerlaw_spr, & !< push-up factor for slip saturation due to twinning constitutive_phenopowerlaw_spr, & !< push-up factor for slip saturation due to twinning
constitutive_phenopowerlaw_twinB, & constitutive_phenopowerlaw_twinB, &
constitutive_phenopowerlaw_twinC, & constitutive_phenopowerlaw_twinC, &
@ -90,19 +85,13 @@ module constitutive_phenopowerlaw
public :: & public :: &
constitutive_phenopowerlaw_init, & constitutive_phenopowerlaw_init, &
#ifndef NEWSTATE
constitutive_phenopowerlaw_stateInit, &
constitutive_phenopowerlaw_aTolState, &
#endif
constitutive_phenopowerlaw_LpAndItsTangent, & constitutive_phenopowerlaw_LpAndItsTangent, &
constitutive_phenopowerlaw_dotState, & constitutive_phenopowerlaw_dotState, &
constitutive_phenopowerlaw_postResults constitutive_phenopowerlaw_postResults
#ifdef NEWSTATE
private :: & private :: &
constitutive_phenopowerlaw_aTolState, & constitutive_phenopowerlaw_aTolState, &
constitutive_phenopowerlaw_stateInit constitutive_phenopowerlaw_stateInit
#endif
contains contains
@ -142,10 +131,8 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
phase_Noutput, & phase_Noutput, &
PLASTICITY_PHENOPOWERLAW_label, & PLASTICITY_PHENOPOWERLAW_label, &
PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_PHENOPOWERLAW_ID, &
material_phase, & material_phase, &
#ifdef NEWSTATE
plasticState, & plasticState, &
#endif
MATERIAL_partPhase MATERIAL_partPhase
use lattice use lattice
use numerics,only: & use numerics,only: &
@ -179,11 +166,8 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifndef NEWSTATE
allocate(constitutive_phenopowerlaw_sizeDotState(maxNinstance), source=0_pInt) allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_sizeState(maxNinstance), source=0_pInt)
#endif
allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_sizePostResult(maxval(phase_Noutput),maxNinstance), & allocate(constitutive_phenopowerlaw_sizePostResult(maxval(phase_Noutput),maxNinstance), &
source=0_pInt) source=0_pInt)
allocate(constitutive_phenopowerlaw_output(maxval(phase_Noutput),maxNinstance)) allocate(constitutive_phenopowerlaw_output(maxval(phase_Noutput),maxNinstance))
@ -428,6 +412,7 @@ allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance),
enddo parsingFile enddo parsingFile
sanityChecks: do phase = 1_pInt, size(phase_plasticity) sanityChecks: do phase = 1_pInt, size(phase_plasticity)
NofMyPhase=count(material_phase==phase)
myPhase: if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then myPhase: if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
constitutive_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance) = & constitutive_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance) = &
@ -486,9 +471,12 @@ allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance),
maxNinstance), source=0.0_pReal) maxNinstance), source=0.0_pReal)
initializeInstances: do phase = 1_pInt, size(phase_plasticity) initializeInstances: do phase = 1_pInt, size(phase_plasticity)
if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then myPhase2: if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
instance = phase_plasticityInstance(phase) instance = phase_plasticityInstance(phase)
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(instance) outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(instance)
select case(constitutive_phenopowerlaw_outputID(o,instance)) select case(constitutive_phenopowerlaw_outputID(o,instance))
case(resistance_slip_ID, & case(resistance_slip_ID, &
@ -515,15 +503,8 @@ allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance),
constitutive_phenopowerlaw_sizePostResults(instance) = constitutive_phenopowerlaw_sizePostResults(instance) + mySize constitutive_phenopowerlaw_sizePostResults(instance) = constitutive_phenopowerlaw_sizePostResults(instance) + mySize
endif outputFound endif outputFound
enddo outputsLoop enddo outputsLoop
!--------------------------------------------------------------------------------------------------
#ifndef NEWSTATE ! allocate state arrays
constitutive_phenopowerlaw_sizeDotState(instance) = constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)+ &
2_pInt + &
constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance) ! s_slip, s_twin, sum(gamma), sum(f), accshear_slip, accshear_twin
constitutive_phenopowerlaw_sizeState(instance) = constitutive_phenopowerlaw_sizeDotState(instance)
#else
sizeState = constitutive_phenopowerlaw_totalNslip(instance)+ & sizeState = constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)+ & constitutive_phenopowerlaw_totalNtwin(instance)+ &
2_pInt + & 2_pInt + &
@ -533,24 +514,23 @@ allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance),
sizeDotState = sizeState sizeDotState = sizeState
plasticState(phase)%sizeDotState = sizeState plasticState(phase)%sizeDotState = sizeState
plasticState(phase)%sizePostResults = constitutive_phenopowerlaw_sizePostResults(instance) plasticState(phase)%sizePostResults = constitutive_phenopowerlaw_sizePostResults(instance)
allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal) allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal)
allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%partionedState0(sizeState,NofMyPhase), source=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)%subState0 (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(plasticState(phase)%dotState_backup(sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)
allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase),source=0.0_pReal)
endif endif
if (any(numerics_integrator == 4_pInt)) & if (any(numerics_integrator == 4_pInt)) &
allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 5_pInt)) & if (any(numerics_integrator == 5_pInt)) &
allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
#endif
do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X
index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1_pInt,instance)) index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1_pInt,instance))
do j = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! loop over (active) systems in my family (slip) do j = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! loop over (active) systems in my family (slip)
@ -601,59 +581,50 @@ allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance),
enddo; enddo enddo; enddo
enddo; enddo enddo; enddo
#ifdef NEWSTATE
call constitutive_phenopowerlaw_stateInit(phase,instance) call constitutive_phenopowerlaw_stateInit(phase,instance)
call constitutive_phenopowerlaw_aTolState(phase,instance) call constitutive_phenopowerlaw_aTolState(phase,instance)
#endif endif myPhase2
endif
enddo initializeInstances enddo initializeInstances
end subroutine constitutive_phenopowerlaw_init end subroutine constitutive_phenopowerlaw_init
#ifdef NEWSTATE
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief sets the initial microstructural state for a given instance of this plasticity !> @brief sets the initial microstructural state for a given instance of this plasticity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_stateInit(phase,instance) subroutine constitutive_phenopowerlaw_stateInit(ph,instance)
use lattice, only: & use lattice, only: &
lattice_maxNslipFamily, & lattice_maxNslipFamily, &
lattice_maxNtwinFamily lattice_maxNtwinFamily
use material, only: & use material, only: &
plasticState plasticState, &
mappingConstitutive
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
instance, & !< number specifying the instance of the plasticity instance, & !< number specifying the instance of the plasticity
phase ph
integer(pInt) :: & integer(pInt) :: &
i i
real(pReal), dimension(size(plasticState(phase)%state(:,1))) :: tempState real(pReal), dimension(plasticState(ph)%sizeState) :: &
tempState
tempState = 0.0_pReal
do i = 1_pInt,lattice_maxNslipFamily do i = 1_pInt,lattice_maxNslipFamily
tempState(1+& tempState(1+sum(constitutive_phenopowerlaw_Nslip(1:i-1,instance)) : &
sum(constitutive_phenopowerlaw_Nslip(1:i-1,instance)) : & sum(constitutive_phenopowerlaw_Nslip(1:i ,instance))) = &
sum(constitutive_phenopowerlaw_Nslip(1:i ,instance))) = & constitutive_phenopowerlaw_tau0_slip(i,instance)
constitutive_phenopowerlaw_tau0_slip(i,instance)
enddo enddo
do i = 1_pInt,lattice_maxNtwinFamily do i = 1_pInt,lattice_maxNtwinFamily
tempState(1+sum(constitutive_phenopowerlaw_Nslip(:,instance))+& tempState(1+sum(constitutive_phenopowerlaw_Nslip(:,instance))+&
sum(constitutive_phenopowerlaw_Ntwin(1:i-1,instance)) : & sum(constitutive_phenopowerlaw_Ntwin(1:i-1,instance)) : &
sum(constitutive_phenopowerlaw_Nslip(:,instance))+& sum(constitutive_phenopowerlaw_Nslip(:,instance))+&
sum(constitutive_phenopowerlaw_Ntwin(1:i ,instance))) = & sum(constitutive_phenopowerlaw_Ntwin(1:i ,instance))) = &
constitutive_phenopowerlaw_tau0_twin(i,instance) constitutive_phenopowerlaw_tau0_twin(i,instance)
enddo enddo
plasticState(phase)%state = spread(tempState,2,size(plasticState(phase)%state(1,:))) plasticState(ph)%state0 = spread(tempState,2,size(plasticState(ph)%state0(1,:)))
plasticState(phase)%state0 = plasticState(phase)%state
plasticState(phase)%partionedState0 = plasticState(phase)%state
end subroutine constitutive_phenopowerlaw_stateInit end subroutine constitutive_phenopowerlaw_stateInit
@ -662,107 +633,37 @@ end subroutine constitutive_phenopowerlaw_stateInit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief sets the relevant state values for a given instance of this plasticity !> @brief sets the relevant state values for a given instance of this plasticity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_aTolState(phase,instance) subroutine constitutive_phenopowerlaw_aTolState(ph,instance)
use material, only: & use material, only: &
plasticState plasticState
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
instance, & !< number specifying the instance of the plasticity instance, & !< number specifying the instance of the plasticity
phase ph
real(pReal), dimension(size(plasticState(phase)%aTolState(:))) :: tempTol
tempTol = 0.0_pReal plasticState(ph)%aTolState(1:constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolResistance(instance)
plasticState(ph)%aTolState(1+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolShear(instance)
plasticState(ph)%aTolState(2+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolTwinFrac(instance)
plasticState(ph)%aTolState(3+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance): &
2+2*(constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance))) = &
constitutive_phenopowerlaw_aTolShear(instance)
tempTol(1:constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolResistance(instance)
tempTol(1+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolShear(instance)
tempTol(2+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolTwinFrac(instance)
tempTol(3+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance): &
2+2*(constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance))) = &
constitutive_phenopowerlaw_aTolShear(instance)
plasticState(phase)%aTolState = tempTol
end subroutine constitutive_phenopowerlaw_aTolState end subroutine constitutive_phenopowerlaw_aTolState
#else
!--------------------------------------------------------------------------------------------------
!> @brief sets the initial microstructural state for a given instance of this plasticity
!--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_stateInit(instance)
use lattice, only: &
lattice_maxNslipFamily, &
lattice_maxNtwinFamily
implicit none
integer(pInt), intent(in) :: &
instance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(instance)) :: &
constitutive_phenopowerlaw_stateInit
integer(pInt) :: &
i
constitutive_phenopowerlaw_stateInit = 0.0_pReal
do i = 1_pInt,lattice_maxNslipFamily
constitutive_phenopowerlaw_stateInit(1+&
sum(constitutive_phenopowerlaw_Nslip(1:i-1,instance)) : &
sum(constitutive_phenopowerlaw_Nslip(1:i ,instance))) = &
constitutive_phenopowerlaw_tau0_slip(i,instance)
enddo
do i = 1_pInt,lattice_maxNtwinFamily
constitutive_phenopowerlaw_stateInit(1+sum(constitutive_phenopowerlaw_Nslip(:,instance))+&
sum(constitutive_phenopowerlaw_Ntwin(1:i-1,instance)) : &
sum(constitutive_phenopowerlaw_Nslip(:,instance))+&
sum(constitutive_phenopowerlaw_Ntwin(1:i ,instance))) = &
constitutive_phenopowerlaw_tau0_twin(i,instance)
enddo
end function constitutive_phenopowerlaw_stateInit
!--------------------------------------------------------------------------------------------------
!> @brief sets the relevant state values for a given instance of this plasticity
!--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_aTolState(instance)
implicit none
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_phenopowerlaw_sizeState(instance)) :: &
constitutive_phenopowerlaw_aTolState
constitutive_phenopowerlaw_aTolState(1:constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolResistance(instance)
constitutive_phenopowerlaw_aTolState(1+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolShear(instance)
constitutive_phenopowerlaw_aTolState(2+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)) = &
constitutive_phenopowerlaw_aTolTwinFrac(instance)
constitutive_phenopowerlaw_aTolState(3+constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance): &
2+2*(constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance))) = &
constitutive_phenopowerlaw_aTolShear(instance)
end function constitutive_phenopowerlaw_aTolState
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,state,ipc,ip,el) subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -784,6 +685,8 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase, & material_phase, &
plasticState, &
mappingConstitutive, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
@ -799,19 +702,13 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
ip, & !< integration point ip, & !< integration point
el !< element el !< element
#ifdef NEWSTATE
real(pReal), dimension(:), intent(in) :: &
state
#else
type(p_vec), intent(in) :: &
state !< microstructure state
#endif
integer(pInt) :: & integer(pInt) :: &
instance, & instance, &
nSlip, & nSlip, &
nTwin,phase,index_Gamma,index_F,index_myFamily, & nTwin,phase,index_Gamma,index_F,index_myFamily, &
f,i,j,k,l,m,n f,i,j,k,l,m,n, &
of, &
ph
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real(pReal), dimension(3,3,2) :: & real(pReal), dimension(3,3,2) :: &
@ -822,8 +719,9 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
gdot_twin,dgdot_dtautwin,tau_twin gdot_twin,dgdot_dtautwin,tau_twin
phase = material_phase(ipc,ip,el) of = mappingConstitutive(1,ipc,ip,el)
instance = phase_plasticityInstance(phase) ph = mappingConstitutive(2,ipc,ip,el)
instance = phase_plasticityInstance(ph)
nSlip = constitutive_phenopowerlaw_totalNslip(instance) nSlip = constitutive_phenopowerlaw_totalNslip(instance)
nTwin = constitutive_phenopowerlaw_totalNtwin(instance) nTwin = constitutive_phenopowerlaw_totalNtwin(instance)
index_Gamma = nSlip + nTwin + 1_pInt index_Gamma = nSlip + nTwin + 1_pInt
@ -835,49 +733,37 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
j = 0_pInt j = 0_pInt
slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j+1_pInt j = j+1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of Lp ! Calculation of Lp
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase)) tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
tau_slip_neg(j) = tau_slip_pos(j) tau_slip_neg(j) = tau_slip_pos(j)
nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,phase) nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1)
do k = 1,lattice_NnonSchmid(phase) do k = 1,lattice_NnonSchmid(ph)
tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* & tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph))
tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* & tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)*& nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,phase) lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)*& nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,phase) lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph)
enddo enddo
#ifdef NEWSTATE
gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* & gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_pos(j))/state(j))**constitutive_phenopowerlaw_n_slip(instance))*& ((abs(tau_slip_pos(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_pos(j)) sign(1.0_pReal,tau_slip_pos(j))
gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* & gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_neg(j))/state(j))**constitutive_phenopowerlaw_n_slip(instance))*& ((abs(tau_slip_neg(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_neg(j)) sign(1.0_pReal,tau_slip_neg(j))
Lp = Lp + (1.0_pReal-state(index_F))*& ! 1-F Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
(gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,phase) (gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph)
#else
gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_pos(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_pos(j))
gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_neg(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_neg(j))
Lp = Lp + (1.0_pReal-state%p(index_F))*& ! 1-F
(gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,phase)
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp ! Calculation of the tangent of Lp
@ -885,7 +771,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
dgdot_dtauslip_pos(j) = gdot_slip_pos(j)*constitutive_phenopowerlaw_n_slip(instance)/tau_slip_pos(j) dgdot_dtauslip_pos(j) = gdot_slip_pos(j)*constitutive_phenopowerlaw_n_slip(instance)/tau_slip_pos(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & 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) + & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip_pos(j)*lattice_Sslip(k,l,1,index_myFamily+i,phase)* & dgdot_dtauslip_pos(j)*lattice_Sslip(k,l,1,index_myFamily+i,ph)* &
nonSchmid_tensor(m,n,1) nonSchmid_tensor(m,n,1)
endif endif
@ -893,7 +779,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
dgdot_dtauslip_neg(j) = gdot_slip_neg(j)*constitutive_phenopowerlaw_n_slip(instance)/tau_slip_neg(j) dgdot_dtauslip_neg(j) = gdot_slip_neg(j)*constitutive_phenopowerlaw_n_slip(instance)/tau_slip_neg(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & 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) + & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip_neg(j)*lattice_Sslip(k,l,1,index_myFamily+i,phase)* & dgdot_dtauslip_neg(j)*lattice_Sslip(k,l,1,index_myFamily+i,ph)* &
nonSchmid_tensor(m,n,2) nonSchmid_tensor(m,n,2)
endif endif
enddo enddo
@ -901,25 +787,18 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
j = 0_pInt j = 0_pInt
twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j+1_pInt j = j+1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of Lp ! Calculation of Lp
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
#ifdef NEWSTATE gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
gdot_twin(j) = (1.0_pReal-state(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*& constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state(nSlip+j))**& (abs(tau_twin(j))/plasticState(ph)%state(nSlip+j,of))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
#else Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,ph)
gdot_twin(j) = (1.0_pReal-state%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
#endif
Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,phase)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp ! Calculation of the tangent of Lp
@ -927,23 +806,21 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(instance)/tau_twin(j) dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(instance)/tau_twin(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & 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) + & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtautwin(j)*lattice_Stwin(k,l,index_myFamily+i,phase)* & dgdot_dtautwin(j)*lattice_Stwin(k,l,index_myFamily+i,ph)* &
lattice_Stwin(m,n,index_myFamily+i,phase) lattice_Stwin(m,n,index_myFamily+i,ph)
endif endif
enddo enddo
enddo twinFamiliesLoop enddo twinFamiliesLoop
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
end subroutine constitutive_phenopowerlaw_LpAndItsTangent
end subroutine constitutive_phenopowerlaw_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) subroutine constitutive_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
use prec, only: &
p_vec
use lattice, only: & use lattice, only: &
lattice_Sslip_v, & lattice_Sslip_v, &
lattice_Stwin_v, & lattice_Stwin_v, &
@ -959,6 +836,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase, & material_phase, &
mappingConstitutive, &
plasticState, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
@ -968,25 +847,14 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element !< microstructure state el !< element !< microstructure state
#ifdef NEWSTATE
real(pReal), dimension(:), intent(in) :: &
state
real(pReal), dimension(size(state)) :: &
constitutive_phenopowerlaw_dotState
#else
type(p_vec), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_dotState
#endif
integer(pInt) :: & integer(pInt) :: &
instance,phase, & instance,ph, &
nSlip,nTwin, & nSlip,nTwin, &
f,i,j,k, & f,i,j,k, &
index_Gamma,index_F,index_myFamily, & index_Gamma,index_F,index_myFamily, &
offset_accshear_slip,offset_accshear_twin offset_accshear_slip,offset_accshear_twin, &
of
real(pReal) :: & real(pReal) :: &
c_SlipSlip,c_SlipTwin,c_TwinSlip,c_TwinTwin, & c_SlipSlip,c_SlipTwin,c_TwinSlip,c_TwinTwin, &
ssat_offset ssat_offset
@ -995,8 +863,10 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
gdot_slip,tau_slip_pos,tau_slip_neg,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip gdot_slip,tau_slip_pos,tau_slip_neg,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,tau_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin gdot_twin,tau_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin
phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase) of = mappingConstitutive(1,ipc,ip,el)
ph = mappingConstitutive(2,ipc,ip,el)
instance = phase_plasticityInstance(ph)
nSlip = constitutive_phenopowerlaw_totalNslip(instance) nSlip = constitutive_phenopowerlaw_totalNslip(instance)
nTwin = constitutive_phenopowerlaw_totalNtwin(instance) nTwin = constitutive_phenopowerlaw_totalNtwin(instance)
@ -1005,89 +875,58 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
index_F = nSlip + nTwin + 2_pInt index_F = nSlip + nTwin + 2_pInt
offset_accshear_slip = nSlip + nTwin + 2_pInt offset_accshear_slip = nSlip + nTwin + 2_pInt
offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip
plasticState(ph)%dotState = 0.0_pReal
constitutive_phenopowerlaw_dotState = 0.0_pReal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
#ifdef NEWSTATE
c_SlipSlip = constitutive_phenopowerlaw_h0_SlipSlip(instance)*& c_SlipSlip = constitutive_phenopowerlaw_h0_SlipSlip(instance)*&
(1.0_pReal + constitutive_phenopowerlaw_twinC(instance)*state(index_F)**& (1.0_pReal + constitutive_phenopowerlaw_twinC(instance)*plasticState(ph)%state(index_F,of)**&
constitutive_phenopowerlaw_twinB(instance)) constitutive_phenopowerlaw_twinB(instance))
c_SlipTwin = 0.0_pReal c_SlipTwin = 0.0_pReal
c_TwinSlip = constitutive_phenopowerlaw_h0_TwinSlip(instance)*& c_TwinSlip = constitutive_phenopowerlaw_h0_TwinSlip(instance)*&
state(index_Gamma)**constitutive_phenopowerlaw_twinE(instance) plasticState(ph)%state(index_Gamma,of)**constitutive_phenopowerlaw_twinE(instance)
c_TwinTwin = constitutive_phenopowerlaw_h0_TwinTwin(instance)*& c_TwinTwin = constitutive_phenopowerlaw_h0_TwinTwin(instance)*&
state(index_F)**constitutive_phenopowerlaw_twinD(instance) plasticState(ph)%state(index_F,of)**constitutive_phenopowerlaw_twinD(instance)
#else
c_SlipSlip = constitutive_phenopowerlaw_h0_SlipSlip(instance)*&
(1.0_pReal + constitutive_phenopowerlaw_twinC(instance)*state%p(index_F)**&
constitutive_phenopowerlaw_twinB(instance))
c_SlipTwin = 0.0_pReal
c_TwinSlip = constitutive_phenopowerlaw_h0_TwinSlip(instance)*&
state%p(index_Gamma)**constitutive_phenopowerlaw_twinE(instance)
c_TwinTwin = constitutive_phenopowerlaw_h0_TwinTwin(instance)*&
state%p(index_F)**constitutive_phenopowerlaw_twinD(instance)
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate left and right vectors and calculate dot gammas ! calculate left and right vectors and calculate dot gammas
#ifdef NEWSTATE ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(plasticState(ph)%state(index_F,of))
ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(state(index_F))
#else
ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(state%p(index_F)) !< microstructure state
#endif
j = 0_pInt j = 0_pInt
slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j+1_pInt j = j+1_pInt
left_SlipSlip(j) = 1.0_pReal ! no system-dependent left part left_SlipSlip(j) = 1.0_pReal ! no system-dependent left part
left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part
#ifdef NEWSTATE right_SlipSlip(j) = abs(1.0_pReal-plasticState(ph)%state(j,of) / &
right_SlipSlip(j) = abs(1.0_pReal-state(j) / &
(constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) & (constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) &
**constitutive_phenopowerlaw_a_slip(instance)& **constitutive_phenopowerlaw_a_slip(instance)&
*sign(1.0_pReal,1.0_pReal-state(j) / & *sign(1.0_pReal,1.0_pReal-plasticState(ph)%state(j,of) / &
(constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) (constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset))
#else
right_SlipSlip(j) = abs(1.0_pReal-state%p(j) / &
(constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) &
**constitutive_phenopowerlaw_a_slip(instance)&
*sign(1.0_pReal,1.0_pReal-state%p(j) / &
(constitutive_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) !< microstructure state
#endif
right_TwinSlip(j) = 1.0_pReal ! no system-dependent part right_TwinSlip(j) = 1.0_pReal ! no system-dependent part
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of dot gamma ! Calculation of dot gamma
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase)) tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
tau_slip_neg(j) = tau_slip_pos(j) tau_slip_neg(j) = tau_slip_pos(j)
do k = 1,lattice_NnonSchmid(phase) do k = 1,lattice_NnonSchmid(ph)
tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* & tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph))
tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* & tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo enddo
#ifdef NEWSTATE
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos(j))/state(j))**constitutive_phenopowerlaw_n_slip(instance) & ((abs(tau_slip_pos(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg(j))/state(j))**constitutive_phenopowerlaw_n_slip(instance))& +(abs(tau_slip_neg(j))/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos(j)) *sign(1.0_pReal,tau_slip_pos(j))
#else
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg(j))/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos(j))
#endif
enddo enddo
enddo slipFamiliesLoop1 enddo slipFamiliesLoop1
j = 0_pInt j = 0_pInt
twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j+1_pInt j = j+1_pInt
left_TwinSlip(j) = 1.0_pReal ! no system-dependent right part left_TwinSlip(j) = 1.0_pReal ! no system-dependent right part
@ -1097,18 +936,11 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Calculation of dot vol frac ! Calculation of dot vol frac
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
#ifdef NEWSTATE gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
gdot_twin(j) = (1.0_pReal-state(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*& constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state(nSlip+j))**& (abs(tau_twin(j))/plasticState(ph)%state(nslip+j,of))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
#else
gdot_twin(j) = (1.0_pReal-state%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
#endif
enddo enddo
enddo twinFamiliesLoop1 enddo twinFamiliesLoop1
@ -1118,54 +950,45 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j+1_pInt j = j+1_pInt
constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j plasticState(ph)%dotState(j,of) = & ! evolution of slip resistance j
c_SlipSlip * left_SlipSlip(j) * & c_SlipSlip * left_SlipSlip(j) * &
dot_product(constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(j,1:nSlip,instance), & dot_product(constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(j,1:nSlip,instance), &
right_SlipSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor right_SlipSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor
c_SlipTwin * left_SlipTwin(j) * & c_SlipTwin * left_SlipTwin(j) * &
dot_product(constitutive_phenopowerlaw_hardeningMatrix_SlipTwin(j,1:nTwin,instance), & dot_product(constitutive_phenopowerlaw_hardeningMatrix_SlipTwin(j,1:nTwin,instance), &
right_SlipTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor right_SlipTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor
constitutive_phenopowerlaw_dotState(index_Gamma) = constitutive_phenopowerlaw_dotState(index_Gamma) + & plasticState(ph)%dotState(index_Gamma,of) = plasticState(ph)%dotState(index_Gamma,of) + &
abs(gdot_slip(j)) abs(gdot_slip(j))
constitutive_phenopowerlaw_dotState(offset_accshear_slip+j) = abs(gdot_slip(j)) plasticState(ph)%dotState(offset_accshear_slip+j,of) = abs(gdot_slip(j))
enddo enddo
enddo slipFamiliesLoop2 enddo slipFamiliesLoop2
j = 0_pInt j = 0_pInt
twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j+1_pInt j = j+1_pInt
constitutive_phenopowerlaw_dotState(j+nSlip) = & ! evolution of twin resistance j plasticState(ph)%dotState(j+nSlip,of) = & ! evolution of twin resistance j
c_TwinSlip * left_TwinSlip(j) * & c_TwinSlip * left_TwinSlip(j) * &
dot_product(constitutive_phenopowerlaw_hardeningMatrix_TwinSlip(j,1:nSlip,instance), & dot_product(constitutive_phenopowerlaw_hardeningMatrix_TwinSlip(j,1:nSlip,instance), &
right_TwinSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor right_TwinSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor
c_TwinTwin * left_TwinTwin(j) * & c_TwinTwin * left_TwinTwin(j) * &
dot_product(constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(j,1:nTwin,instance), & dot_product(constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(j,1:nTwin,instance), &
right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor
#ifndef NEWSTATE if (plasticState(ph)%state(index_F,of) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0
if (state%p(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 plasticState(ph)%dotState(index_F,of) = plasticState(ph)%dotState(index_F,of) + &
constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + & gdot_twin(j)/lattice_shearTwin(index_myFamily+i,ph)
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,phase) plasticState(ph)%dotState(offset_accshear_twin+j,of) = abs(gdot_twin(j))
#else
if (state(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0
constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + &
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,phase)
#endif
constitutive_phenopowerlaw_dotState(offset_accshear_twin+j) = abs(gdot_twin(j))
enddo enddo
enddo twinFamiliesLoop2 enddo twinFamiliesLoop2
end function constitutive_phenopowerlaw_dotState end subroutine constitutive_phenopowerlaw_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el) function constitutive_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -1174,6 +997,8 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains, &
material_phase, & material_phase, &
plasticState, &
mappingConstitutive, &
phase_plasticityInstance, & phase_plasticityInstance, &
phase_Noutput phase_Noutput
use lattice, only: & use lattice, only: &
@ -1195,27 +1020,21 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element !< microstructure state el !< element !< microstructure state
#ifdef NEWSTATE
real(pReal), dimension(:), intent(in) :: &
state
#else
type(p_vec), intent(in) :: &
state !< microstructure state
#endif
real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_postResults constitutive_phenopowerlaw_postResults
integer(pInt) :: & integer(pInt) :: &
instance,phase, & instance,ph, of, &
nSlip,nTwin, & nSlip,nTwin, &
o,f,i,c,j,k, & o,f,i,c,j,k, &
index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily
real(pReal) :: & real(pReal) :: &
tau_slip_pos,tau_slip_neg,tau tau_slip_pos,tau_slip_neg,tau
phase = material_phase(ipc,ip,el) of = mappingConstitutive(1,ipc,ip,el)
instance = phase_plasticityInstance(phase) ph = mappingConstitutive(2,ipc,ip,el)
instance = phase_plasticityInstance(ph)
nSlip = constitutive_phenopowerlaw_totalNslip(instance) nSlip = constitutive_phenopowerlaw_totalNslip(instance)
nTwin = constitutive_phenopowerlaw_totalNtwin(instance) nTwin = constitutive_phenopowerlaw_totalNtwin(instance)
@ -1231,48 +1050,33 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(instance) outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(instance)
select case(constitutive_phenopowerlaw_outputID(o,instance)) select case(constitutive_phenopowerlaw_outputID(o,instance))
case (resistance_slip_ID) case (resistance_slip_ID)
#ifdef NEWSTATE constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = plasticState(ph)%state(1:nSlip,of)
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(1:nSlip)
#else
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state%p(1:nSlip)
#endif
c = c + nSlip c = c + nSlip
case (accumulatedshear_slip_ID) case (accumulatedshear_slip_ID)
#ifdef NEWSTATE constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = plasticState(ph)%state(index_accshear_slip:&
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(index_accshear_slip:& index_accshear_slip+nSlip-1_pInt,of)
index_accshear_slip+nSlip-1_pInt)
#else
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state%p(index_accshear_slip:&
index_accshear_slip+nSlip-1_pInt)
#endif
c = c + nSlip c = c + nSlip
case (shearrate_slip_ID) case (shearrate_slip_ID)
j = 0_pInt j = 0_pInt
slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j + 1_pInt j = j + 1_pInt
tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase)) tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
tau_slip_neg = tau_slip_pos tau_slip_neg = tau_slip_pos
do k = 1,lattice_NnonSchmid(phase) do k = 1,lattice_NnonSchmid(ph)
tau_slip_pos = tau_slip_pos + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* & tau_slip_pos = tau_slip_pos + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph))
tau_slip_neg = tau_slip_neg + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* & tau_slip_neg = tau_slip_neg + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph))
enddo enddo
#ifdef NEWSTATE
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/state(j))**constitutive_phenopowerlaw_n_slip(instance) & ((abs(tau_slip_pos)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/state(j))**constitutive_phenopowerlaw_n_slip(instance))& +(abs(tau_slip_neg)/plasticState(ph)%state(j,of))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos) *sign(1.0_pReal,tau_slip_pos)
#else
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/state%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
+(abs(tau_slip_neg)/state%p(j))**constitutive_phenopowerlaw_n_slip(instance))&
*sign(1.0_pReal,tau_slip_pos)
#endif
enddo enddo
enddo slipFamiliesLoop1 enddo slipFamiliesLoop1
c = c + nSlip c = c + nSlip
@ -1280,63 +1084,40 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
case (resolvedstress_slip_ID) case (resolvedstress_slip_ID)
j = 0_pInt j = 0_pInt
slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j + 1_pInt j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = & constitutive_phenopowerlaw_postResults(c+j) = &
dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
enddo enddo
enddo slipFamiliesLoop2 enddo slipFamiliesLoop2
c = c + nSlip c = c + nSlip
case (totalshear_ID) case (totalshear_ID)
#ifdef NEWSTATE
constitutive_phenopowerlaw_postResults(c+1_pInt) = & constitutive_phenopowerlaw_postResults(c+1_pInt) = &
state(index_Gamma) plasticState(ph)%state(index_Gamma,of)
#else
constitutive_phenopowerlaw_postResults(c+1_pInt) = &
state%p(index_Gamma)
#endif
c = c + 1_pInt c = c + 1_pInt
case (resistance_twin_ID) case (resistance_twin_ID)
#ifdef NEWSTATE
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = & constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
state(1_pInt+nSlip:nTwin+nSlip-1_pInt) plasticState(ph)%state(1_pInt+nSlip:nTwin+nSlip-1_pInt,of)
#else
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
state%p(1_pInt+nSlip:nTwin+nSlip-1_pInt)
#endif
c = c + nTwin c = c + nTwin
case (accumulatedshear_twin_ID) case (accumulatedshear_twin_ID)
#ifdef NEWSTATE
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = & constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
state(index_accshear_twin:index_accshear_twin+nTwin-1_pInt) plasticState(ph)%state(index_accshear_twin:index_accshear_twin+nTwin-1_pInt,of)
#else
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
state%p(index_accshear_twin:index_accshear_twin+nTwin-1_pInt)
#endif
c = c + nTwin c = c + nTwin
case (shearrate_twin_ID) case (shearrate_twin_ID)
j = 0_pInt j = 0_pInt
twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j + 1_pInt j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase)) tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
#ifdef NEWSTATE constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*& constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau)/state(j+nSlip))**& (abs(tau)/plasticState(ph)%state(j+nSlip,of))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau)) constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau))
#else
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau)/state%p(j+nSlip))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau))
#endif
enddo enddo
enddo twinFamiliesLoop1 enddo twinFamiliesLoop1
c = c + nTwin c = c + nTwin
@ -1344,21 +1125,17 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
case (resolvedstress_twin_ID) case (resolvedstress_twin_ID)
j = 0_pInt j = 0_pInt
twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j + 1_pInt j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = & constitutive_phenopowerlaw_postResults(c+j) = &
dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase)) dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph))
enddo enddo
enddo twinFamiliesLoop2 enddo twinFamiliesLoop2
c = c + nTwin c = c + nTwin
case (totalvolfrac_ID) case (totalvolfrac_ID)
#ifdef NEWSTATE constitutive_phenopowerlaw_postResults(c+1_pInt) = plasticState(ph)%state(index_F,of)
constitutive_phenopowerlaw_postResults(c+1_pInt) = state(index_F)
#else
constitutive_phenopowerlaw_postResults(c+1_pInt) = state%p(index_F)
#endif
c = c + 1_pInt c = c + 1_pInt
end select end select
@ -1366,6 +1143,4 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
end function constitutive_phenopowerlaw_postResults end function constitutive_phenopowerlaw_postResults
end module constitutive_phenopowerlaw end module constitutive_phenopowerlaw

23
code/constitutive_thermal.f90
View File

@ -19,7 +19,6 @@ module constitutive_thermal
constitutive_thermal_init, & constitutive_thermal_init, &
constitutive_thermal_microstructure, & constitutive_thermal_microstructure, &
constitutive_thermal_collectDotState, & constitutive_thermal_collectDotState, &
constitutive_thermal_collectDeltaState, &
constitutive_thermal_postResults constitutive_thermal_postResults
contains contains
@ -204,28 +203,6 @@ subroutine constitutive_thermal_collectDotState(Tstar_v, Lp, ipc, ip, el)
end subroutine constitutive_thermal_collectDotState end subroutine constitutive_thermal_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state (so far, only nonlocal)
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
logical function constitutive_thermal_collectDeltaState(ipc, ip, el)
use material, only: &
material_phase, &
phase_thermal
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
select case (phase_thermal(material_phase(ipc,ip,el)))
end select
end function constitutive_thermal_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

570
code/constitutive_titanmod.f90
View File

File diff suppressed because it is too large Load Diff

692
code/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

10
code/damage_gradient.f90
View File

@ -14,8 +14,6 @@ module damage_gradient
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
damage_gradient_sizeDotState, & !< number of dotStates
damage_gradient_sizeState, & !< total number of microstructural state variables
damage_gradient_sizePostResults !< cumulative size of post results damage_gradient_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -112,8 +110,6 @@ subroutine damage_gradient_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(damage_gradient_sizeDotState(maxNinstance), source=0_pInt)
allocate(damage_gradient_sizeState(maxNinstance), source=0_pInt)
allocate(damage_gradient_sizePostResults(maxNinstance), source=0_pInt) allocate(damage_gradient_sizePostResults(maxNinstance), source=0_pInt)
allocate(damage_gradient_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) allocate(damage_gradient_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(damage_gradient_output(maxval(phase_Noutput),maxNinstance)) allocate(damage_gradient_output(maxval(phase_Noutput),maxNinstance))
@ -168,8 +164,6 @@ subroutine damage_gradient_init(fileUnit)
if (phase_damage(phase) == DAMAGE_gradient_ID) then if (phase_damage(phase) == DAMAGE_gradient_ID) then
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
instance = phase_damageInstance(phase) instance = phase_damageInstance(phase)
damage_gradient_sizeDotState(instance) = 1_pInt
damage_gradient_sizeState(instance) = 3_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Determine size of postResults array ! Determine size of postResults array
@ -187,8 +181,8 @@ subroutine damage_gradient_init(fileUnit)
endif endif
enddo outputsLoop enddo outputsLoop
! Determine size of state array ! Determine size of state array
sizeDotState = damage_gradient_sizeDotState(instance) sizeDotState = 1_pInt
sizeState = damage_gradient_sizeState (instance) sizeState = 3_pInt
damageState(phase)%sizeState = sizeState damageState(phase)%sizeState = sizeState
damageState(phase)%sizeDotState = sizeDotState damageState(phase)%sizeDotState = sizeDotState

10
code/damage_none.f90
View File

@ -12,8 +12,6 @@ module damage_none
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
damage_none_sizeDotState, &
damage_none_sizeState, &
damage_none_sizePostResults damage_none_sizePostResults
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -69,10 +67,9 @@ subroutine damage_none_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifdef NEWSTATE
initializeInstances: do phase = 1_pInt, size(phase_damage) initializeInstances: do phase = 1_pInt, size(phase_damage)
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
if (phase_damage(phase) == DAMAGE_none_ID .and. NofMyPhase/=0) then if (phase_damage(phase) == DAMAGE_none_ID) then
sizeState = 0_pInt sizeState = 0_pInt
damageState(phase)%sizeState = sizeState damageState(phase)%sizeState = sizeState
sizeDotState = sizeState sizeDotState = sizeState
@ -95,10 +92,7 @@ subroutine damage_none_init(fileUnit)
allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase)) allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase))
endif endif
enddo initializeInstances enddo initializeInstances
#else
allocate(damage_none_sizeDotState(maxNinstance), source=1_pInt)
allocate(damage_none_sizeState(maxNinstance), source=1_pInt)
#endif
allocate(damage_none_sizePostResults(maxNinstance), source=0_pInt) allocate(damage_none_sizePostResults(maxNinstance), source=0_pInt)
end subroutine damage_none_init end subroutine damage_none_init

36
code/homogenization.f90
View File

@ -104,12 +104,10 @@ subroutine homogenization_init()
FE_geomtype FE_geomtype
use constitutive, only: & use constitutive, only: &
constitutive_maxSizePostResults constitutive_maxSizePostResults
#ifdef NEWSTATE
use constitutive_damage, only: & use constitutive_damage, only: &
constitutive_damage_maxSizePostResults constitutive_damage_maxSizePostResults
use constitutive_thermal, only: & use constitutive_thermal, only: &
constitutive_thermal_maxSizePostResults constitutive_thermal_maxSizePostResults
#endif
use crystallite, only: & use crystallite, only: &
crystallite_maxSizePostResults crystallite_maxSizePostResults
use material use material
@ -239,10 +237,8 @@ subroutine homogenization_init()
materialpoint_sizeResults = 1 & ! grain count materialpoint_sizeResults = 1 & ! grain count
+ 1 + homogenization_maxSizePostResults & ! homogSize & homogResult + 1 + homogenization_maxSizePostResults & ! homogSize & homogResult
+ homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results + homogenization_maxNgrains * (1 + crystallite_maxSizePostResults & ! crystallite size & crystallite results
#ifdef NEWSTATE
+ 1 + constitutive_damage_maxSizePostResults & + 1 + constitutive_damage_maxSizePostResults &
+ 1 + constitutive_thermal_maxSizePostResults & + 1 + constitutive_thermal_maxSizePostResults &
#endif
+ 1 + constitutive_maxSizePostResults) ! constitutive size & constitutive results + 1 + constitutive_maxSizePostResults) ! constitutive size & constitutive results
allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems)) allocate(materialpoint_results(materialpoint_sizeResults,mesh_maxNips,mesh_NcpElems))
@ -302,19 +298,12 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
#ifdef NEWSTATE
plasticState, & plasticState, &
damageState, & damageState, &
thermalState, & thermalState, &
mappingConstitutive, & mappingConstitutive, &
#endif
homogenization_Ngrains homogenization_Ngrains
#ifndef NEWSTATE
use constitutive, only: &
constitutive_state0, &
constitutive_partionedState0, &
constitutive_state
#endif
use crystallite, only: & use crystallite, only: &
crystallite_heat, & crystallite_heat, &
@ -380,16 +369,12 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains) forall(i = FEsolving_execIP(1,e):FEsolving_execIP(2,e), g = 1:myNgrains)
#ifdef NEWSTATE
plasticState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = & plasticState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%state0(:,mappingConstitutive(1,g,i,e)) plasticState(mappingConstitutive(2,g,i,e))%state0(:,mappingConstitutive(1,g,i,e))
damageState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = & damageState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
damageState(mappingConstitutive(2,g,i,e))%state0(:,mappingConstitutive(1,g,i,e)) damageState(mappingConstitutive(2,g,i,e))%state0(:,mappingConstitutive(1,g,i,e))
thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = & thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%state0(:,mappingConstitutive(1,g,i,e)) thermalState(mappingConstitutive(2,g,i,e))%state0(:,mappingConstitutive(1,g,i,e))
#else
constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures
#endif
crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads crystallite_partionedFp0(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) ! ...plastic def grads
crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads crystallite_partionedLp0(1:3,1:3,g,i,e) = crystallite_Lp0(1:3,1:3,g,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness
@ -444,7 +429,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e)! ...stiffness crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e)! ...stiffness
crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) = crystallite_Tstar_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
#ifdef NEWSTATE
forall (g = 1:myNgrains) forall (g = 1:myNgrains)
plasticState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = & plasticState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e))
@ -453,10 +437,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = & thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) thermalState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e))
end forall end forall
#else
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructures
#endif
if (homogenization_sizeState(i,e) > 0_pInt) & if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme homogenization_subState0(i,e)%p = homogenization_state(i,e)%p ! ...internal state of homog scheme
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e) ! ...def grad
@ -503,7 +483,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness crystallite_dPdF(1:3,1:3,1:3,1:3,1:myNgrains,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,1:myNgrains,i,e) ! ...stiffness
crystallite_Tstar_v(1:6,1:myNgrains,i,e) = crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress crystallite_Tstar_v(1:6,1:myNgrains,i,e) = crystallite_partionedTstar0_v(1:6,1:myNgrains,i,e) ! ...2nd PK stress
#ifdef NEWSTATE
forall (g = 1:myNgrains) forall (g = 1:myNgrains)
plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = & plasticState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = &
plasticState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) plasticState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e))
@ -512,10 +491,6 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
thermalState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = & thermalState(mappingConstitutive(2,g,i,e))%state(:,mappingConstitutive(1,g,i,e)) = &
thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e)) thermalState(mappingConstitutive(2,g,i,e))%partionedState0(:,mappingConstitutive(1,g,i,e))
end forall end forall
#else
forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures
#endif
if (homogenization_sizeState(i,e) > 0_pInt) & if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_state(i,e)%p = homogenization_subState0(i,e)%p ! ...internal state of homog scheme homogenization_state(i,e)%p = homogenization_subState0(i,e)%p ! ...internal state of homog scheme
endif endif
@ -631,18 +606,13 @@ subroutine materialpoint_postResults
use mesh, only: & use mesh, only: &
mesh_element mesh_element
use material, only: & use material, only: &
#ifdef NEWSTATE
plasticState, & plasticState, &
damageState, & damageState, &
thermalState, & thermalState, &
material_phase, & material_phase, &
#endif
homogenization_Ngrains, & homogenization_Ngrains, &
microstructure_crystallite microstructure_crystallite
use constitutive, only: & use constitutive, only: &
#ifndef NEWSTATE
constitutive_sizePostResults, &
#endif
constitutive_postResults constitutive_postResults
use crystallite, only: & use crystallite, only: &
crystallite_sizePostResults, & crystallite_sizePostResults, &
@ -679,13 +649,9 @@ subroutine materialpoint_postResults
grainLooping :do g = 1,myNgrains grainLooping :do g = 1,myNgrains
theSize = (1 + crystallite_sizePostResults(myCrystallite)) + & theSize = (1 + crystallite_sizePostResults(myCrystallite)) + &
#ifdef NEWSTATE
(1 + plasticState(material_phase(g,i,e))%sizePostResults) + & (1 + plasticState(material_phase(g,i,e))%sizePostResults) + &
(1 + damageState(material_phase(g,i,e))%sizePostResults) + & (1 + damageState(material_phase(g,i,e))%sizePostResults) + &
(1 + thermalState(material_phase(g,i,e))%sizePostResults) (1 + thermalState(material_phase(g,i,e))%sizePostResults)
#else
(1 + constitutive_sizePostResults(g,i,e))
#endif
materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e) ! tell crystallite results materialpoint_results(thePos+1:thePos+theSize,i,e) = crystallite_postResults(g,i,e) ! tell crystallite results
thePos = thePos + theSize thePos = thePos + theSize
enddo grainLooping enddo grainLooping

11
code/lattice.f90
View File

@ -635,14 +635,12 @@ module lattice
real(pReal), dimension(:), allocatable, public, protected :: & real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, & lattice_mu, &
lattice_nu lattice_nu
#ifdef NEWSTATE
real(pReal), dimension(:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_thermalConductivity33, & lattice_thermalConductivity33, &
lattice_thermalExpansion33, & lattice_thermalExpansion33, &
lattice_surfaceEnergy33 lattice_surfaceEnergy33
real(pReal), dimension(:), allocatable, public, protected :: & real(pReal), dimension(:), allocatable, public, protected :: &
lattice_referenceTemperature lattice_referenceTemperature
#endif
enum, bind(c) enum, bind(c)
enumerator :: LATTICE_undefined_ID, & enumerator :: LATTICE_undefined_ID, &
LATTICE_iso_ID, & LATTICE_iso_ID, &
@ -864,12 +862,10 @@ subroutine lattice_init
allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID) allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID)
allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) allocate(lattice_C66(6,6,Nphases), source=0.0_pReal)
allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal) allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal)
#ifdef NEWSTATE
allocate(lattice_thermalConductivity33(3,3,Nphases), source=0.0_pReal) allocate(lattice_thermalConductivity33(3,3,Nphases), source=0.0_pReal)
allocate(lattice_thermalExpansion33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_thermalExpansion33 (3,3,Nphases), source=0.0_pReal)
allocate(lattice_surfaceEnergy33 (3,3,Nphases), source=0.0_pReal) allocate(lattice_surfaceEnergy33 (3,3,Nphases), source=0.0_pReal)
allocate(lattice_referenceTemperature (Nphases), source=0.0_pReal) allocate(lattice_referenceTemperature (Nphases), source=0.0_pReal)
#endif
allocate(lattice_mu(Nphases), source=0.0_pReal) allocate(lattice_mu(Nphases), source=0.0_pReal)
allocate(lattice_nu(Nphases), source=0.0_pReal) allocate(lattice_nu(Nphases), source=0.0_pReal)
@ -955,7 +951,6 @@ subroutine lattice_init
lattice_C66(6,6,section) = IO_floatValue(line,positions,2_pInt) lattice_C66(6,6,section) = IO_floatValue(line,positions,2_pInt)
case ('covera_ratio','c/a_ratio','c/a') case ('covera_ratio','c/a_ratio','c/a')
CoverA(section) = IO_floatValue(line,positions,2_pInt) CoverA(section) = IO_floatValue(line,positions,2_pInt)
#ifdef NEWSTATE
case ('k11') case ('k11')
lattice_thermalConductivity33(1,1,section) = IO_floatValue(line,positions,2_pInt) lattice_thermalConductivity33(1,1,section) = IO_floatValue(line,positions,2_pInt)
case ('k22') case ('k22')
@ -976,7 +971,6 @@ subroutine lattice_init
lattice_surfaceEnergy33(3,3,section) = IO_floatValue(line,positions,2_pInt) lattice_surfaceEnergy33(3,3,section) = IO_floatValue(line,positions,2_pInt)
case ('reference_temperature') case ('reference_temperature')
lattice_referenceTemperature(section) = IO_floatValue(line,positions,2_pInt) lattice_referenceTemperature(section) = IO_floatValue(line,positions,2_pInt)
#endif
end select end select
endif endif
enddo enddo
@ -1048,14 +1042,12 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
do i = 1_pInt, 6_pInt do i = 1_pInt, 6_pInt
if (abs(lattice_C66(i,i,myPhase))<tol_math_check) call IO_error(43_pInt,el=i,ip=myPhase) if (abs(lattice_C66(i,i,myPhase))<tol_math_check) call IO_error(43_pInt,el=i,ip=myPhase)
enddo enddo
#ifdef NEWSTATE
lattice_thermalConductivity33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_thermalConductivity33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_thermalConductivity33(1:3,1:3,myPhase)) lattice_thermalConductivity33(1:3,1:3,myPhase))
lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_thermalExpansion33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_thermalExpansion33(1:3,1:3,myPhase)) lattice_thermalExpansion33(1:3,1:3,myPhase))
lattice_surfaceEnergy33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),& lattice_surfaceEnergy33(1:3,1:3,myPhase) = lattice_symmetrize33(lattice_structure(myPhase),&
lattice_surfaceEnergy33(1:3,1:3,myPhase)) lattice_surfaceEnergy33(1:3,1:3,myPhase))
#endif
select case(lattice_structure(myPhase)) select case(lattice_structure(myPhase))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1267,8 +1259,6 @@ pure function lattice_symmetrizeC66(struct,C66)
end function lattice_symmetrizeC66 end function lattice_symmetrizeC66
#ifdef NEWSTATE
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Symmetrizes 2nd order tensor according to lattice type !> @brief Symmetrizes 2nd order tensor according to lattice type
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1298,7 +1288,6 @@ pure function lattice_symmetrize33(struct,T33)
end select end select
end function lattice_symmetrize33 end function lattice_symmetrize33
#endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

30
code/material.f90
View File

@ -12,9 +12,7 @@ module material
use prec, only: & use prec, only: &
pReal, & pReal, &
pInt, & pInt, &
#ifdef NEWSTATE
tState, & tState, &
#endif
p_intvec p_intvec
implicit none implicit none
@ -51,7 +49,6 @@ module material
PLASTICITY_titanmod_ID, & PLASTICITY_titanmod_ID, &
PLASTICITY_nonlocal_ID PLASTICITY_nonlocal_ID
end enum end enum
#ifdef NEWSTATE
enum, bind(c) enum, bind(c)
enumerator :: DAMAGE_none_ID, & enumerator :: DAMAGE_none_ID, &
DAMAGE_local_ID, & DAMAGE_local_ID, &
@ -63,7 +60,6 @@ module material
THERMAL_conduction_ID, & THERMAL_conduction_ID, &
THERMAL_adiabatic_ID THERMAL_adiabatic_ID
end enum end enum
#endif
enum, bind(c) enum, bind(c)
enumerator :: HOMOGENIZATION_undefined_ID, & enumerator :: HOMOGENIZATION_undefined_ID, &
HOMOGENIZATION_none_ID, & HOMOGENIZATION_none_ID, &
@ -84,12 +80,11 @@ module material
phase_elasticity !< elasticity of each phase phase_elasticity !< elasticity of each phase
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: &
phase_plasticity !< plasticity of each phase phase_plasticity !< plasticity of each phase
#ifdef NEWSTATE
integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: &
phase_damage !< damage of each phase phase_damage !< damage of each phase
integer(kind(THERMAL_none_ID)), dimension(:), allocatable, public, protected :: & integer(kind(THERMAL_none_ID)), dimension(:), allocatable, public, protected :: &
phase_thermal !< thermal of each phase phase_thermal !< thermal of each phase
#endif
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: &
homogenization_type !< type of each homogenization homogenization_type !< type of each homogenization
@ -111,10 +106,8 @@ module material
phase_Noutput, & !< number of '(output)' items per phase phase_Noutput, & !< number of '(output)' items per phase
phase_elasticityInstance, & !< instance of particular elasticity of each phase phase_elasticityInstance, & !< instance of particular elasticity of each phase
phase_plasticityInstance, & !< instance of particular plasticity of each phase phase_plasticityInstance, & !< instance of particular plasticity of each phase
#ifdef NEWSTATE phase_damageInstance, & !< instance of particular damage of each phase
phase_damageInstance, & !< instance of particular plasticity of each phase phase_thermalInstance, & !< instance of particular thermal of each phase
phase_thermalInstance, & !< instance of particular plasticity of each phase
#endif
crystallite_Noutput, & !< number of '(output)' items per crystallite setting crystallite_Noutput, & !< number of '(output)' items per crystallite setting
homogenization_typeInstance, & !< instance of particular type of each homogenization homogenization_typeInstance, & !< instance of particular type of each homogenization
microstructure_crystallite !< crystallite setting ID of each microstructure microstructure_crystallite !< crystallite setting ID of each microstructure
@ -122,13 +115,12 @@ module material
integer(pInt), dimension(:,:,:), allocatable, public :: & integer(pInt), dimension(:,:,:), allocatable, public :: &
material_phase !< phase (index) of each grain,IP,element material_phase !< phase (index) of each grain,IP,element
#ifdef NEWSTATE
type(tState), allocatable, dimension(:), public :: & type(tState), allocatable, dimension(:), public :: &
plasticState, & plasticState, &
elasticState, & elasticState, &
damageState, & damageState, &
thermalState thermalState
#endif
integer(pInt), dimension(:,:,:), allocatable, public, protected :: & integer(pInt), dimension(:,:,:), allocatable, public, protected :: &
material_texture !< texture (index) of each grain,IP,element material_texture !< texture (index) of each grain,IP,element
@ -181,12 +173,11 @@ module material
logical, dimension(:), allocatable, private :: & logical, dimension(:), allocatable, private :: &
homogenization_active homogenization_active
#if defined(HDF) || defined(NEWSTATE)
integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: mappingConstitutive integer(pInt), dimension(:,:,:,:), allocatable, public, protected :: mappingConstitutive
integer(pInt), dimension(:,:,:), allocatable, public, protected :: mappingCrystallite integer(pInt), dimension(:,:,:), allocatable, public, protected :: mappingCrystallite
integer(pInt), dimension(:), allocatable :: ConstitutivePosition integer(pInt), dimension(:), allocatable :: ConstitutivePosition
integer(pInt), dimension(:), allocatable :: CrystallitePosition integer(pInt), dimension(:), allocatable :: CrystallitePosition
#endif
public :: & public :: &
material_init, & material_init, &
@ -197,7 +188,6 @@ module material
PLASTICITY_dislotwin_ID, & PLASTICITY_dislotwin_ID, &
PLASTICITY_titanmod_ID, & PLASTICITY_titanmod_ID, &
PLASTICITY_nonlocal_ID, & PLASTICITY_nonlocal_ID, &
#ifdef NEWSTATE
DAMAGE_none_ID, & DAMAGE_none_ID, &
DAMAGE_local_ID, & DAMAGE_local_ID, &
DAMAGE_gradient_ID, & DAMAGE_gradient_ID, &
@ -205,7 +195,6 @@ module material
THERMAL_iso_ID, & THERMAL_iso_ID, &
THERMAL_conduction_ID, & THERMAL_conduction_ID, &
THERMAL_adiabatic_ID, & THERMAL_adiabatic_ID, &
#endif
HOMOGENIZATION_none_ID, & HOMOGENIZATION_none_ID, &
HOMOGENIZATION_isostrain_ID, & HOMOGENIZATION_isostrain_ID, &
#ifdef HDF #ifdef HDF
@ -276,12 +265,11 @@ subroutine material_init
call material_parsePhase(FILEUNIT,material_partPhase) call material_parsePhase(FILEUNIT,material_partPhase)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed' if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'
close(FILEUNIT) close(FILEUNIT)
#ifdef NEWSTATE
allocate(plasticState(material_Nphase)) allocate(plasticState(material_Nphase))
allocate(elasticState(material_Nphase)) allocate(elasticState(material_Nphase))
allocate(damageState (material_Nphase)) allocate(damageState (material_Nphase))
allocate(thermalState(material_Nphase)) allocate(thermalState(material_Nphase))
#endif
do m = 1_pInt,material_Nmicrostructure do m = 1_pInt,material_Nmicrostructure
if(microstructure_crystallite(m) < 1_pInt .or. & if(microstructure_crystallite(m) < 1_pInt .or. &
@ -321,7 +309,6 @@ subroutine material_init
call material_populateGrains call material_populateGrains
#if defined(HDF) || defined(NEWSTATE)
allocate(mappingConstitutive(2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt) allocate(mappingConstitutive(2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
allocate(mappingCrystallite (2,homogenization_maxNgrains,mesh_NcpElems),source=0_pInt) allocate(mappingCrystallite (2,homogenization_maxNgrains,mesh_NcpElems),source=0_pInt)
allocate(ConstitutivePosition(material_Nphase),source=0_pInt) allocate(ConstitutivePosition(material_Nphase),source=0_pInt)
@ -335,7 +322,6 @@ subroutine material_init
enddo GrainLoop enddo GrainLoop
enddo IPloop enddo IPloop
enddo ElemLoop enddo ElemLoop
#endif
end subroutine material_init end subroutine material_init
@ -632,12 +618,10 @@ subroutine material_parsePhase(fileUnit,myPart)
allocate(phase_elasticityInstance(Nsections), source=0_pInt) allocate(phase_elasticityInstance(Nsections), source=0_pInt)
allocate(phase_plasticity(Nsections) , source=PLASTICITY_undefined_ID) allocate(phase_plasticity(Nsections) , source=PLASTICITY_undefined_ID)
allocate(phase_plasticityInstance(Nsections), source=0_pInt) allocate(phase_plasticityInstance(Nsections), source=0_pInt)
#ifdef NEWSTATE
allocate(phase_damage(Nsections) , source=DAMAGE_none_ID) allocate(phase_damage(Nsections) , source=DAMAGE_none_ID)
allocate(phase_damageInstance(Nsections), source=0_pInt) allocate(phase_damageInstance(Nsections), source=0_pInt)
allocate(phase_thermal(Nsections) , source=THERMAL_none_ID) allocate(phase_thermal(Nsections) , source=THERMAL_none_ID)
allocate(phase_thermalInstance(Nsections), source=0_pInt) allocate(phase_thermalInstance(Nsections), source=0_pInt)
#endif
allocate(phase_Noutput(Nsections), source=0_pInt) allocate(phase_Noutput(Nsections), source=0_pInt)
allocate(phase_localPlasticity(Nsections), source=.false.) allocate(phase_localPlasticity(Nsections), source=.false.)
@ -694,7 +678,6 @@ subroutine material_parsePhase(fileUnit,myPart)
call IO_error(201_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt))) call IO_error(201_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt)))
end select end select
phase_plasticityInstance(section) = count(phase_plasticity == phase_plasticity(section)) ! count instances phase_plasticityInstance(section) = count(phase_plasticity == phase_plasticity(section)) ! count instances
#ifdef NEWSTATE
case ('damage') case ('damage')
select case (IO_lc(IO_stringValue(line,positions,2_pInt))) select case (IO_lc(IO_stringValue(line,positions,2_pInt)))
case (DAMAGE_NONE_label) case (DAMAGE_NONE_label)
@ -721,7 +704,6 @@ subroutine material_parsePhase(fileUnit,myPart)
call IO_error(200_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt))) call IO_error(200_pInt,ext_msg=trim(IO_stringValue(line,positions,2_pInt)))
end select end select
phase_thermalInstance(section) = count(phase_thermal == phase_thermal(section)) ! count instances phase_thermalInstance(section) = count(phase_thermal == phase_thermal(section)) ! count instances
#endif
end select end select
endif endif
enddo enddo

9
code/prec.f90
View File

@ -57,10 +57,11 @@ module prec
integer(pInt), dimension(:), allocatable :: p integer(pInt), dimension(:), allocatable :: p
end type p_intvec end type p_intvec
#ifdef NEWSTATE
!http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array !http://stackoverflow.com/questions/3948210/can-i-have-a-pointer-to-an-item-in-an-allocatable-array
type, public :: tState type, public :: tState
integer(pInt) :: sizeState,sizeDotState,sizePostResults integer(pInt) :: sizeState = 0_pInt , &
sizeDotState = 0_pInt, &
sizePostResults = 0_pInt
logical :: nonlocal logical :: nonlocal
real(pReal), allocatable, dimension(:) :: atolState real(pReal), allocatable, dimension(:) :: atolState
real(pReal), allocatable, dimension(:,:) :: state, & ! material points, state size real(pReal), allocatable, dimension(:,:) :: state, & ! material points, state size
@ -76,7 +77,6 @@ module prec
RK4dotState RK4dotState
real(pReal), allocatable, dimension(:,:,:) :: RKCK45dotState real(pReal), allocatable, dimension(:,:,:) :: RKCK45dotState
end type end type
#endif
public :: & public :: &
prec_init prec_init
@ -96,9 +96,6 @@ subroutine prec_init
write(6,'(/,a)') ' <<<+- prec init -+>>>' write(6,'(/,a)') ' <<<+- prec init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
#include "compilation_info.f90" #include "compilation_info.f90"
#ifdef NEWSTATE
write(6,'(a)') 'Using new state structure'
#endif
write(6,'(a,i3)') ' Bytes for pReal: ',pReal write(6,'(a,i3)') ' Bytes for pReal: ',pReal
write(6,'(a,i3)') ' Bytes for pInt: ',pInt write(6,'(a,i3)') ' Bytes for pInt: ',pInt
write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt write(6,'(a,i3)') ' Bytes for pLongInt: ',pLongInt

11
code/thermal_adiabatic.f90
View File

@ -14,8 +14,6 @@ module thermal_adiabatic
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
thermal_adiabatic_sizeDotState, & !< number of dotStates
thermal_adiabatic_sizeState, & !< total number of microstructural state variables
thermal_adiabatic_sizePostResults !< cumulative size of post results thermal_adiabatic_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -111,8 +109,6 @@ subroutine thermal_adiabatic_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(thermal_adiabatic_sizeDotState(maxNinstance), source=0_pInt)
allocate(thermal_adiabatic_sizeState(maxNinstance), source=0_pInt)
allocate(thermal_adiabatic_sizePostResults(maxNinstance), source=0_pInt) allocate(thermal_adiabatic_sizePostResults(maxNinstance), source=0_pInt)
allocate(thermal_adiabatic_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) allocate(thermal_adiabatic_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(thermal_adiabatic_output(maxval(phase_Noutput),maxNinstance)) allocate(thermal_adiabatic_output(maxval(phase_Noutput),maxNinstance))
@ -165,8 +161,6 @@ subroutine thermal_adiabatic_init(fileUnit)
if (phase_thermal(phase) == THERMAL_adiabatic_ID) then if (phase_thermal(phase) == THERMAL_adiabatic_ID) then
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
instance = phase_thermalInstance(phase) instance = phase_thermalInstance(phase)
thermal_adiabatic_sizeDotState(instance) = 1_pInt
thermal_adiabatic_sizeState(instance) = 1_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Determine size of postResults array ! Determine size of postResults array
@ -182,9 +176,8 @@ subroutine thermal_adiabatic_init(fileUnit)
endif endif
enddo outputsLoop enddo outputsLoop
! Determine size of state array ! Determine size of state array
sizeDotState = thermal_adiabatic_sizeDotState(instance) sizeDotState = 1_pInt
sizeState = thermal_adiabatic_sizeState (instance) sizeState = 1_pInt
thermalState(phase)%sizeState = sizeState thermalState(phase)%sizeState = sizeState
thermalState(phase)%sizeDotState = sizeDotState thermalState(phase)%sizeDotState = sizeDotState
allocate(thermalState(phase)%aTolState (sizeState), source=0.0_pReal) allocate(thermalState(phase)%aTolState (sizeState), source=0.0_pReal)

12
code/thermal_conduction.f90
View File

@ -14,8 +14,6 @@ module thermal_conduction
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
thermal_conduction_sizeDotState, & !< number of dotStates
thermal_conduction_sizeState, & !< total number of microstructural state variables
thermal_conduction_sizePostResults !< cumulative size of post results thermal_conduction_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -111,8 +109,6 @@ subroutine thermal_conduction_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
allocate(thermal_conduction_sizeDotState(maxNinstance), source=0_pInt)
allocate(thermal_conduction_sizeState(maxNinstance), source=0_pInt)
allocate(thermal_conduction_sizePostResults(maxNinstance), source=0_pInt) allocate(thermal_conduction_sizePostResults(maxNinstance), source=0_pInt)
allocate(thermal_conduction_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) allocate(thermal_conduction_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
allocate(thermal_conduction_output(maxval(phase_Noutput),maxNinstance)) allocate(thermal_conduction_output(maxval(phase_Noutput),maxNinstance))
@ -165,8 +161,6 @@ subroutine thermal_conduction_init(fileUnit)
if (phase_thermal(phase) == THERMAL_conduction_ID) then if (phase_thermal(phase) == THERMAL_conduction_ID) then
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
instance = phase_thermalInstance(phase) instance = phase_thermalInstance(phase)
thermal_conduction_sizeDotState(instance) = 0_pInt
thermal_conduction_sizeState(instance) = 2_pInt
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Determine size of postResults array ! Determine size of postResults array
@ -182,8 +176,8 @@ subroutine thermal_conduction_init(fileUnit)
endif endif
enddo outputsLoop enddo outputsLoop
! Determine size of state array ! Determine size of state array
sizeDotState = thermal_conduction_sizeDotState(instance) sizeDotState = 0_pInt
sizeState = thermal_conduction_sizeState (instance) sizeState = 2_pInt
thermalState(phase)%sizeState = sizeState thermalState(phase)%sizeState = sizeState
thermalState(phase)%sizeDotState = sizeDotState thermalState(phase)%sizeDotState = sizeDotState
@ -195,7 +189,7 @@ subroutine thermal_conduction_init(fileUnit)
allocate(thermalState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) allocate(thermalState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal)
allocate(thermalState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(thermalState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(thermalState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(thermalState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal)
allocate(thermalState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(thermalState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal)
if (any(numerics_integrator == 1_pInt)) then if (any(numerics_integrator == 1_pInt)) then
allocate(thermalState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) allocate(thermalState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal)

9
code/thermal_none.f90
View File

@ -12,8 +12,6 @@ module thermal_none
implicit none implicit none
private private
integer(pInt), dimension(:), allocatable, public, protected :: & integer(pInt), dimension(:), allocatable, public, protected :: &
thermal_none_sizeDotState, &
thermal_none_sizeState, &
thermal_none_sizePostResults thermal_none_sizePostResults
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
@ -68,8 +66,7 @@ subroutine thermal_none_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
#ifdef NEWSTATE
initializeInstances: do phase = 1_pInt, size(phase_thermal) initializeInstances: do phase = 1_pInt, size(phase_thermal)
NofMyPhase=count(material_phase==phase) NofMyPhase=count(material_phase==phase)
if (phase_thermal(phase) == THERMAL_none_ID .and. NofMyPhase/=0) then if (phase_thermal(phase) == THERMAL_none_ID .and. NofMyPhase/=0) then
@ -95,10 +92,6 @@ subroutine thermal_none_init(fileUnit)
allocate(thermalState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase)) allocate(thermalState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase))
endif endif
enddo initializeInstances enddo initializeInstances
#else
allocate(thermal_none_sizeDotState(maxNinstance), source=1_pInt)
allocate(thermal_none_sizeState(maxNinstance), source=1_pInt)
#endif
allocate(thermal_none_sizePostResults(maxNinstance), source=0_pInt) allocate(thermal_none_sizePostResults(maxNinstance), source=0_pInt)
end subroutine thermal_none_init end subroutine thermal_none_init