Update in order to include the c/a-ratio as additional information about the structure (modules lattice, constitutive_ changed accordingly)

Correction of the expression for the velocity of mobile dislocation densities in constitutive_dislobased.f90

Correction of some Fortran syntax, incompatible with Compaq Visual Fortran (e.g. integer(8))
This commit is contained in:
Luc Hantcherli 2009-03-20 14:34:24 +00:00
parent 121723ffa5
commit e5255f4384
11 changed files with 758 additions and 624 deletions

View File

@ -485,17 +485,24 @@
use prec, only: pReal,pInt,reltol_Outer use prec, only: pReal,pInt,reltol_Outer
use constitutive, only: constitutive_dotState,constitutive_sizeDotState,& use constitutive, only: constitutive_dotState,constitutive_sizeDotState,&
constitutive_state_old,constitutive_state_new constitutive_state_old,constitutive_state_new
use debug
logical integrateState logical integrateState
integer(pInt) tick,tock,tickrate,maxticks
integer(pInt) g,i,e,mySize integer(pInt) g,i,e,mySize
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6) :: Tstar_v
real(pReal) dt real(pReal) dt
real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum
mySize = constitutive_sizeDotState(g,i,e) mySize = constitutive_sizeDotState(g,i,e)
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
residuum = constitutive_state_new(g,i,e)%p(1:mySize) - constitutive_state_old(g,i,e)%p(1:mySize) - & residuum = constitutive_state_new(g,i,e)%p(1:mySize) - constitutive_state_old(g,i,e)%p(1:mySize) - &
dt*constitutive_dotState(Tstar_v,CPFEM_Temperature(i,e),g,i,e) ! residuum from evolution of microstructure dt*constitutive_dotState(Tstar_v,CPFEM_Temperature(i,e),g,i,e) ! residuum from evolution of microstructure
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick
if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks
constitutive_state_new(g,i,e)%p(1:mySize) = constitutive_state_new(g,i,e)%p(1:mySize) - residuum ! update of microstructure constitutive_state_new(g,i,e)%p(1:mySize) = constitutive_state_new(g,i,e)%p(1:mySize) - residuum ! update of microstructure
integrateState = maxval(abs(residuum/constitutive_state_new(g,i,e)%p(1:mySize)),& integrateState = maxval(abs(residuum/constitutive_state_new(g,i,e)%p(1:mySize)),&
constitutive_state_new(g,i,e)%p(1:mySize) /= 0.0_pReal) < reltol_Outer constitutive_state_new(g,i,e)%p(1:mySize) /= 0.0_pReal) < reltol_Outer
@ -652,7 +659,7 @@
logical failed logical failed
integer(pInt) cp_en, ip, grain integer(pInt) cp_en, ip, grain
integer(pInt) iInner,dummy, i,j,k,l,m,n integer(pInt) iInner,dummy, i,j,k,l,m,n
integer(8) tick,tock,tickrate,maxticks integer(pInt) tick,tock,tickrate,maxticks
real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2 real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
@ -722,7 +729,7 @@ Inner: do ! inner iteration: Lp
! check for acceleration/deceleration in Newton--Raphson correction ! check for acceleration/deceleration in Newton--Raphson correction
! !
if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed
leapfrog == 1.0) then leapfrog == 1.0) then
Lpguess = Lpguess_old ! restore known good guess Lpguess = Lpguess_old ! restore known good guess
msg = 'NaN present' ! croak for cutback msg = 'NaN present' ! croak for cutback
return return

View File

@ -40,7 +40,6 @@
character(256) path character(256) path
inquire(6, name=path) ! determine outputfile inquire(6, name=path) ! determine outputfile
write(6,*)'trying to open mattex at',path(1:scan(path,pathSep,back=.true.))//relPath
open(unit,status='old',err=100,file=path(1:scan(path,pathSep,back=.true.))//relPath) open(unit,status='old',err=100,file=path(1:scan(path,pathSep,back=.true.))//relPath)
IO_open_file = .true. IO_open_file = .true.
return return

View File

@ -45,6 +45,7 @@ MODULE constitutive_dislobased
constitutive_dislobased_sizeState, & constitutive_dislobased_sizeState, &
constitutive_dislobased_sizePostResults constitutive_dislobased_sizePostResults
character(len=64), dimension(:,:), allocatable :: constitutive_dislobased_output character(len=64), dimension(:,:), allocatable :: constitutive_dislobased_output
character(len=32), dimension(:), allocatable :: constitutive_dislobased_structureName
integer(pInt), dimension(:), allocatable :: constitutive_dislobased_structure integer(pInt), dimension(:), allocatable :: constitutive_dislobased_structure
integer(pInt), dimension(:), allocatable :: constitutive_dislobased_Nslip integer(pInt), dimension(:), allocatable :: constitutive_dislobased_Nslip
real(pReal), dimension(:), allocatable :: constitutive_dislobased_C11 real(pReal), dimension(:), allocatable :: constitutive_dislobased_C11
@ -68,6 +69,7 @@ MODULE constitutive_dislobased
real(pReal), dimension(:), allocatable :: constitutive_dislobased_c6 real(pReal), dimension(:), allocatable :: constitutive_dislobased_c6
real(pReal), dimension(:), allocatable :: constitutive_dislobased_c7 real(pReal), dimension(:), allocatable :: constitutive_dislobased_c7
real(pReal), dimension(:), allocatable :: constitutive_dislobased_c8 real(pReal), dimension(:), allocatable :: constitutive_dislobased_c8
real(pReal), dimension(:), allocatable :: constitutive_dislobased_CoverA
real(pReal), dimension(:,:), allocatable :: constitutive_dislobased_SlipIntCoeff real(pReal), dimension(:,:), allocatable :: constitutive_dislobased_SlipIntCoeff
real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Iparallel real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Iparallel
real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Iforest real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Iforest
@ -102,7 +104,7 @@ subroutine constitutive_dislobased_init(file)
use math, only: math_Mandel3333to66, math_Voigt66to3333, math_mul3x3 use math, only: math_Mandel3333to66, math_Voigt66to3333, math_mul3x3
use IO use IO
use material use material
use lattice, only: lattice_sn, lattice_st, lattice_SlipIntType use lattice, only: lattice_sn, lattice_st, lattice_interactionSlipSlip, lattice_initializeStructure
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7 integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
@ -120,6 +122,7 @@ subroutine constitutive_dislobased_init(file)
allocate(constitutive_dislobased_sizePostResults(maxNinstance)); constitutive_dislobased_sizePostResults = 0_pInt allocate(constitutive_dislobased_sizePostResults(maxNinstance)); constitutive_dislobased_sizePostResults = 0_pInt
allocate(constitutive_dislobased_output(maxval(phase_Noutput), & allocate(constitutive_dislobased_output(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_dislobased_output = '' maxNinstance)) ; constitutive_dislobased_output = ''
allocate(constitutive_dislobased_structureName(maxNinstance)) ; constitutive_dislobased_structureName = ''
allocate(constitutive_dislobased_structure(maxNinstance)) ; constitutive_dislobased_structure = 0_pInt allocate(constitutive_dislobased_structure(maxNinstance)) ; constitutive_dislobased_structure = 0_pInt
allocate(constitutive_dislobased_Nslip(maxNinstance)) ; constitutive_dislobased_Nslip = 0_pInt allocate(constitutive_dislobased_Nslip(maxNinstance)) ; constitutive_dislobased_Nslip = 0_pInt
allocate(constitutive_dislobased_C11(maxNinstance)) ; constitutive_dislobased_C11 = 0.0_pReal allocate(constitutive_dislobased_C11(maxNinstance)) ; constitutive_dislobased_C11 = 0.0_pReal
@ -142,6 +145,7 @@ subroutine constitutive_dislobased_init(file)
allocate(constitutive_dislobased_c6(maxNinstance)) ; constitutive_dislobased_c6 = 0.0_pReal allocate(constitutive_dislobased_c6(maxNinstance)) ; constitutive_dislobased_c6 = 0.0_pReal
allocate(constitutive_dislobased_c7(maxNinstance)) ; constitutive_dislobased_c7 = 0.0_pReal allocate(constitutive_dislobased_c7(maxNinstance)) ; constitutive_dislobased_c7 = 0.0_pReal
allocate(constitutive_dislobased_c8(maxNinstance)) ; constitutive_dislobased_c8 = 0.0_pReal allocate(constitutive_dislobased_c8(maxNinstance)) ; constitutive_dislobased_c8 = 0.0_pReal
allocate(constitutive_dislobased_CoverA(maxNinstance)) ; constitutive_dislobased_CoverA = 0.0_pReal
allocate(constitutive_dislobased_SlipIntCoeff(6,maxNinstance)) ; constitutive_dislobased_SlipIntCoeff = 0.0_pReal allocate(constitutive_dislobased_SlipIntCoeff(6,maxNinstance)) ; constitutive_dislobased_SlipIntCoeff = 0.0_pReal
rewind(file) rewind(file)
@ -169,7 +173,9 @@ subroutine constitutive_dislobased_init(file)
output = output + 1 output = output + 1
constitutive_dislobased_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) constitutive_dislobased_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('lattice_structure') case ('lattice_structure')
constitutive_dislobased_structure(i) = IO_intValue(line,positions,2) constitutive_dislobased_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
case ('covera_ratio')
constitutive_dislobased_CoverA(i) = IO_floatValue(line,positions,2)
case ('nslip') case ('nslip')
constitutive_dislobased_Nslip(i) = IO_intValue(line,positions,2) constitutive_dislobased_Nslip(i) = IO_intValue(line,positions,2)
case ('c11') case ('c11')
@ -208,17 +214,19 @@ subroutine constitutive_dislobased_init(file)
constitutive_dislobased_c7(i) = IO_floatValue(line,positions,2) constitutive_dislobased_c7(i) = IO_floatValue(line,positions,2)
case ('c8') case ('c8')
constitutive_dislobased_c8(i) = IO_floatValue(line,positions,2) constitutive_dislobased_c8(i) = IO_floatValue(line,positions,2)
case ('interaction_coefficients') case ('interaction_coefficients')
forall (j=2:min(7,positions(1))) & forall (j=2:min(7,positions(1))) &
constitutive_dislobased_SlipIntCoeff(j-1,i) = IO_floatValue(line,positions,j) constitutive_dislobased_SlipIntCoeff(j-1,i) = IO_floatValue(line,positions,j)
end select end select
endif endif
enddo enddo
100 do i = 1,maxNinstance ! sanity checks 100 do i = 1,maxNinstance
if (constitutive_dislobased_structure(i) < 1 .or. & constitutive_dislobased_structure(i) = lattice_initializeStructure(constitutive_dislobased_structureName(i), &
constitutive_dislobased_structure(i) > 3) call IO_error(201) constitutive_dislobased_CoverA(i))
! sanity checks
if (constitutive_dislobased_structure(i) < 1) call IO_error(201)
if (constitutive_dislobased_Nslip(i) < 1) call IO_error(202) if (constitutive_dislobased_Nslip(i) < 1) call IO_error(202)
if (constitutive_dislobased_rho0(i) < 0.0_pReal) call IO_error(220) if (constitutive_dislobased_rho0(i) < 0.0_pReal) call IO_error(220)
if (constitutive_dislobased_bg(i) <= 0.0_pReal) call IO_error(221) if (constitutive_dislobased_bg(i) <= 0.0_pReal) call IO_error(221)
@ -259,7 +267,7 @@ subroutine constitutive_dislobased_init(file)
constitutive_dislobased_Cslip_66(k,k,i) = constitutive_dislobased_C11(i) constitutive_dislobased_Cslip_66(k,k,i) = constitutive_dislobased_C11(i)
constitutive_dislobased_Cslip_66(k+3,k+3,i) = constitutive_dislobased_C44(i) constitutive_dislobased_Cslip_66(k+3,k+3,i) = constitutive_dislobased_C44(i)
end forall end forall
case(3) ! hcp case(3:) ! all hex
constitutive_dislobased_Cslip_66(1,1,i) = constitutive_dislobased_C11(i) constitutive_dislobased_Cslip_66(1,1,i) = constitutive_dislobased_C11(i)
constitutive_dislobased_Cslip_66(2,2,i) = constitutive_dislobased_C11(i) constitutive_dislobased_Cslip_66(2,2,i) = constitutive_dislobased_C11(i)
constitutive_dislobased_Cslip_66(3,3,i) = constitutive_dislobased_C33(i) constitutive_dislobased_Cslip_66(3,3,i) = constitutive_dislobased_C33(i)
@ -286,14 +294,11 @@ subroutine constitutive_dislobased_init(file)
x = math_mul3x3(lattice_sn(:,j,i),lattice_st(:,k,i)) x = math_mul3x3(lattice_sn(:,j,i),lattice_st(:,k,i))
y = 1.0_pReal-x**(2.0_pReal) y = 1.0_pReal-x**(2.0_pReal)
!* Interaction matrix * !* Interaction matrix *
constitutive_dislobased_Iforest(j,k,i)=abs(x)*& constitutive_dislobased_Iforest(j,k,i) = abs(x)*&
constitutive_dislobased_SlipIntCoeff(lattice_SlipIntType(j,k,constitutive_dislobased_structure(i)),i) constitutive_dislobased_SlipIntCoeff(lattice_interactionSlipSlip(j,k,constitutive_dislobased_structure(i)),i)
if (y>0.0_pReal) then if (y>0.0_pReal) &
constitutive_dislobased_Iparallel(j,k,i)=sqrt(y)*& constitutive_dislobased_Iparallel(j,k,i) = sqrt(y)*&
constitutive_dislobased_SlipIntCoeff(lattice_SlipIntType(j,k,constitutive_dislobased_structure(i)),i) constitutive_dislobased_SlipIntCoeff(lattice_interactionSlipSlip(j,k,constitutive_dislobased_structure(i)),i)
else
constitutive_dislobased_Iparallel(j,k,i)=0.0_pReal
endif
enddo enddo
enddo enddo
@ -303,7 +308,7 @@ subroutine constitutive_dislobased_init(file)
end subroutine end subroutine
function constitutive_dislobased_stateInit(ipc,ip,el) function constitutive_dislobased_stateInit(ipc,ip,el)
!********************************************************************* !*********************************************************************
!* initial microstructural state * !* initial microstructural state *
@ -407,7 +412,7 @@ subroutine constitutive_dislobased_microstructure(Temperature,state,ipc,ip,el)
state(ipc,ip,el)%p(7*n+i) = & state(ipc,ip,el)%p(7*n+i) = &
state(ipc,ip,el)%p(6*n+i)*constitutive_dislobased_bg(matID)*attack_frequency*state(ipc,ip,el)%p(4*n+i)*& state(ipc,ip,el)%p(6*n+i)*constitutive_dislobased_bg(matID)*attack_frequency*state(ipc,ip,el)%p(4*n+i)*&
exp(-constitutive_dislobased_Qedge(matID)/(kB*Temperature)) exp(-constitutive_dislobased_Qedge(matID)/(kB*Temperature))
enddo enddo
end subroutine end subroutine
@ -450,20 +455,24 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera
!* Calculation of Lp !* Calculation of Lp
Lp = 0.0_pReal Lp = 0.0_pReal
gdot_slip = 0.0_pReal
do i = 1,constitutive_dislobased_Nslip(matID) do i = 1,constitutive_dislobased_Nslip(matID)
tau_slip(i) = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) tau_slip(i) = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID)))
gdot_slip(i) = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip(i))*& if ((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))>0) &
sinh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) gdot_slip(i) = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip(i))*&
sinh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature))
Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_dislobased_structure(matID)) Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_dislobased_structure(matID))
enddo enddo
!* Calculation of the tangent of Lp !* Calculation of the tangent of Lp
dLp_dTstar3333 = 0.0_pReal dLp_dTstar3333 = 0.0_pReal
dLp_dTstar = 0.0_pReal dLp_dTstar = 0.0_pReal
dgdot_dtauslip = 0.0_pReal
do i = 1,constitutive_dislobased_Nslip(matID) do i = 1,constitutive_dislobased_Nslip(matID)
dgdot_dtauslip(i) = (state(ipc,ip,el)%p(7*n+i)*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)*& if ((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))>0) &
cosh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) dgdot_dtauslip(i) = (state(ipc,ip,el)%p(7*n+i)*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)*&
cosh(((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature))
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
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(i)*lattice_Sslip(k,l,i,constitutive_dislobased_structure(matID))* & dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_dislobased_structure(matID))* &
@ -511,13 +520,17 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el)
if (abs(tau_slip) > state(ipc,ip,el)%p(3*n+i)) then if (abs(tau_slip) > state(ipc,ip,el)%p(3*n+i)) then
gdot_slip = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip)*& gdot_slip = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip)*&
sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature))
locks = (sqrt(state(ipc,ip,el)%p(n+i))*abs(gdot_slip))/& locks = (sqrt(state(ipc,ip,el)%p(n+i))*abs(gdot_slip))/&
(constitutive_dislobased_c4(matID)*constitutive_dislobased_bg(matID)) (constitutive_dislobased_c4(matID)*constitutive_dislobased_bg(matID))
athermal_recovery = constitutive_dislobased_c7(matID)*state(ipc,ip,el)%p(i)*abs(gdot_slip) athermal_recovery = constitutive_dislobased_c7(matID)*state(ipc,ip,el)%p(i)*abs(gdot_slip)
thermal_recovery = constitutive_dislobased_c8(matID)*abs(tau_slip)*state(ipc,ip,el)%p(i)**(2.0_pReal)*&
((constitutive_dislobased_D0(matID)*constitutive_dislobased_bg(matID)**(3.0_pReal))/& !thermal_recovery = constitutive_dislobased_c8(matID)*abs(tau_slip)*state(ipc,ip,el)%p(i)**(2.0_pReal)*&
(kB*Temperature))*exp(-constitutive_dislobased_Qsd(matID)/(kB*Temperature)) ! ((constitutive_dislobased_D0(matID)*constitutive_dislobased_bg(matID)**(3.0_pReal))/&
constitutive_dislobased_dotState(i) = locks - athermal_recovery !-thermal_recovery ! (kB*Temperature))*exp(-constitutive_dislobased_Qsd(matID)/(kB*Temperature))
constitutive_dislobased_dotState(i) = locks - athermal_recovery
endif endif
enddo enddo
@ -564,8 +577,12 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i
case ('rateofshear') case ('rateofshear')
do i = 1,n do i = 1,n
tau_slip = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) tau_slip = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID)))
constitutive_dislobased_postResults(c+i) = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip)*& if ((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))>0) then
sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature)) constitutive_dislobased_postResults(c+i) = state(ipc,ip,el)%p(7*n+i)*sign(1.0_pReal,tau_slip)*&
sinh(((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))*state(ipc,ip,el)%p(5*n+i))/(kB*Temperature))
else
constitutive_dislobased_postResults(c+i) = 0.0_pReal
endif
enddo enddo
c = c + n c = c + n
end select end select

View File

@ -298,18 +298,18 @@ subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,sta
if (norm_Tstar > 0) then if (norm_Tstar > 0) then
!* Calculation of Lp !* Calculation of Lp
Lp = Tstar33/norm_Tstar*constitutive_j2_gdot0(matID)/constitutive_j2_fTaylor(matID)* & Lp = Tstar33/norm_Tstar*constitutive_j2_gdot0(matID)/constitutive_j2_fTaylor(matID)* &
(dsqrt(1.5_pReal)/constitutive_j2_fTaylor(matID)*norm_Tstar/state(ipc,ip,el)%p(1))**constitutive_j2_n(matID) (dsqrt(1.5_pReal)/constitutive_j2_fTaylor(matID)*norm_Tstar/state(ipc,ip,el)%p(1))**constitutive_j2_n(matID)
!* Calculation of the tangent of Lp !* Calculation of the tangent of Lp
factor = constitutive_j2_gdot0(matID)/constitutive_j2_fTaylor(matID)* & factor = constitutive_j2_gdot0(matID)/constitutive_j2_fTaylor(matID)* &
(dsqrt(1.5_pReal)/ constitutive_j2_fTaylor(matID)/state(ipc,ip,el)%p(1))**constitutive_j2_n(matID) * & (dsqrt(1.5_pReal)/ constitutive_j2_fTaylor(matID)/state(ipc,ip,el)%p(1))**constitutive_j2_n(matID) * &
norm_Tstar**(constitutive_j2_n(matID)-1.0_pReal) norm_Tstar**(constitutive_j2_n(matID)-1.0_pReal)
dLp_dTstar3333 = 0.0_pReal dLp_dTstar3333 = 0.0_pReal
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dTstar3333(k,l,m,n) = Tstar33(k,l)*Tstar33(m,n) * (constitutive_j2_n(matID)-1.0_pReal)/squarenorm_Tstar dLp_dTstar3333(k,l,m,n) = Tstar33(k,l)*Tstar33(m,n) * (constitutive_j2_n(matID)-1.0_pReal)/squarenorm_Tstar
forall (k=1:3,l=1:3) & forall (k=1:3,l=1:3) &
dLp_dTstar3333(k,l,k,l) = dLp_dTstar3333(k,l,k,l) + 1.0_pReal dLp_dTstar3333(k,l,k,l) = dLp_dTstar3333(k,l,k,l) + 1.0_pReal
dLp_dTstar = math_Plain3333to99(factor * dLp_dTstar3333) dLp_dTstar = math_Plain3333to99(factor * dLp_dTstar3333)
end if end if
return return

View File

@ -37,9 +37,11 @@ MODULE constitutive_phenomenological
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_sizeDotState, & integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_sizeDotState, &
constitutive_phenomenological_sizeState, & constitutive_phenomenological_sizeState, &
constitutive_phenomenological_sizePostResults constitutive_phenomenological_sizePostResults
character(len=64), dimension(:,:), allocatable :: constitutive_phenomenological_output character(len=64), dimension(:,:), allocatable :: constitutive_phenomenological_output
character(len=32), dimension(:), allocatable :: constitutive_phenomenological_structureName
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_structure integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_structure
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_Nslip integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_Nslip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_CoverA
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C11 real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C11
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C12 real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C12
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C13 real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C13
@ -77,7 +79,8 @@ subroutine constitutive_phenomenological_init(file)
use prec, only: pInt, pReal use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333 use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO use IO
use material use material
use lattice, only: lattice_initializeStructure
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7 integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
@ -92,9 +95,11 @@ subroutine constitutive_phenomenological_init(file)
allocate(constitutive_phenomenological_sizeState(maxNinstance)) ; constitutive_phenomenological_sizeState = 0_pInt allocate(constitutive_phenomenological_sizeState(maxNinstance)) ; constitutive_phenomenological_sizeState = 0_pInt
allocate(constitutive_phenomenological_sizePostResults(maxNinstance)); constitutive_phenomenological_sizePostResults = 0_pInt allocate(constitutive_phenomenological_sizePostResults(maxNinstance)); constitutive_phenomenological_sizePostResults = 0_pInt
allocate(constitutive_phenomenological_output(maxval(phase_Noutput), & allocate(constitutive_phenomenological_output(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_phenomenological_output = '' maxNinstance)) ; constitutive_phenomenological_output = ''
allocate(constitutive_phenomenological_structureName(maxNinstance)) ; constitutive_phenomenological_structureName = ''
allocate(constitutive_phenomenological_structure(maxNinstance)) ; constitutive_phenomenological_structure = 0_pInt allocate(constitutive_phenomenological_structure(maxNinstance)) ; constitutive_phenomenological_structure = 0_pInt
allocate(constitutive_phenomenological_Nslip(maxNinstance)) ; constitutive_phenomenological_Nslip = 0_pInt allocate(constitutive_phenomenological_Nslip(maxNinstance)) ; constitutive_phenomenological_Nslip = 0_pInt
allocate(constitutive_phenomenological_CoverA(maxNinstance)) ; constitutive_phenomenological_CoverA = 0.0_pReal
allocate(constitutive_phenomenological_C11(maxNinstance)) ; constitutive_phenomenological_C11 = 0.0_pReal allocate(constitutive_phenomenological_C11(maxNinstance)) ; constitutive_phenomenological_C11 = 0.0_pReal
allocate(constitutive_phenomenological_C12(maxNinstance)) ; constitutive_phenomenological_C12 = 0.0_pReal allocate(constitutive_phenomenological_C12(maxNinstance)) ; constitutive_phenomenological_C12 = 0.0_pReal
allocate(constitutive_phenomenological_C13(maxNinstance)) ; constitutive_phenomenological_C13 = 0.0_pReal allocate(constitutive_phenomenological_C13(maxNinstance)) ; constitutive_phenomenological_C13 = 0.0_pReal
@ -133,10 +138,12 @@ subroutine constitutive_phenomenological_init(file)
case ('(output)') case ('(output)')
output = output + 1 output = output + 1
constitutive_phenomenological_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) constitutive_phenomenological_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('lattice_structure') case ('lattice_structure')
constitutive_phenomenological_structure(i) = IO_intValue(line,positions,2) constitutive_phenomenological_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
case ('nslip') case ('nslip')
constitutive_phenomenological_Nslip(i) = IO_intValue(line,positions,2) constitutive_phenomenological_Nslip(i) = IO_intValue(line,positions,2)
case ('covera_ratio')
constitutive_phenomenological_CoverA(i) = IO_floatValue(line,positions,2)
case ('c11') case ('c11')
constitutive_phenomenological_C11(i) = IO_floatValue(line,positions,2) constitutive_phenomenological_C11(i) = IO_floatValue(line,positions,2)
case ('c12') case ('c12')
@ -165,7 +172,9 @@ subroutine constitutive_phenomenological_init(file)
endif endif
enddo enddo
100 do i = 1,maxNinstance ! sanity checks 100 do i = 1,maxNinstance
constitutive_phenomenological_structure(i) = lattice_initializeStructure(constitutive_phenomenological_structureName(i), &
constitutive_phenomenological_CoverA(i)) ! sanity checks
if (constitutive_phenomenological_structure(i) < 1 .or. & if (constitutive_phenomenological_structure(i) < 1 .or. &
constitutive_phenomenological_structure(i) > 3) call IO_error(201) constitutive_phenomenological_structure(i) > 3) call IO_error(201)
if (constitutive_phenomenological_Nslip(i) < 1) call IO_error(202) if (constitutive_phenomenological_Nslip(i) < 1) call IO_error(202)

View File

@ -9,8 +9,10 @@
integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution = 0_pInt integer(pInt), dimension(nCutback+1) :: debug_cutbackDistribution = 0_pInt
integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution = 0_pInt integer(pInt), dimension(nInner) :: debug_InnerLoopDistribution = 0_pInt
integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution = 0_pInt integer(pInt), dimension(nOuter) :: debug_OuterLoopDistribution = 0_pInt
integer(8) :: debug_cumLpTicks = 0_pInt integer(pInt) :: debug_cumLpTicks = 0_pInt
integer(pInt) :: debug_cumDotStateTicks = 0_pInt
integer(pInt) :: debug_cumLpCalls = 0_pInt integer(pInt) :: debug_cumLpCalls = 0_pInt
integer(pInt) :: debug_cumDotStateCalls = 0_pInt
logical :: debugger = .false. logical :: debugger = .false.
logical :: distribution_init = .false. logical :: distribution_init = .false.
@ -26,18 +28,26 @@
implicit none implicit none
integer(pInt) i,integral integer(pInt) i,integral
integer(8) tickrate integer(pInt) tickrate
write(6,*) write(6,*)
write(6,*) 'DEBUG Info' write(6,*) 'DEBUG Info'
write(6,*) write(6,*)
write(6,'(a33,x,i9)') 'total calls to LpAndItsTangent :',debug_cumLpCalls write(6,'(a33,x,i12)') 'total calls to LpAndItsTangent :',debug_cumLpCalls
if (debug_cumLpCalls > 0_pInt) then if (debug_cumLpCalls > 0_pInt) then
call system_clock(count_rate=tickrate) call system_clock(count_rate=tickrate)
write(6,'(a33,x,f10.6)') 'avg CPU time/microsecs per call :',dble(debug_cumLpTicks)/tickrate/1.0e-6_pReal/debug_cumLpCalls write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',dble(debug_cumLpTicks)/tickrate/1.0e-6_pReal/debug_cumLpCalls
write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumLpTicks write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumLpTicks
endif endif
write(6,*) write(6,*)
write(6,'(a33,x,i12)') 'total calls to dotState :',debug_cumDotStateCalls
if (debug_cumdotStateCalls > 0_pInt) then
call system_clock(count_rate=tickrate)
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotStateTicks)/tickrate/1.0e-6_pReal/debug_cumDotStateCalls
write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks
endif
write(6,*)
write(6,*) 'distribution_cutback :' write(6,*) 'distribution_cutback :'
do i=0,nCutback do i=0,nCutback
if (debug_cutbackDistribution(i+1) /= 0) write(6,*) i,debug_cutbackDistribution(i+1) if (debug_cutbackDistribution(i+1) /= 0) write(6,*) i,debug_cutbackDistribution(i+1)

File diff suppressed because it is too large Load Diff

View File

@ -13,6 +13,10 @@ use prec, only: pReal,pInt
implicit none implicit none
character(len=64), parameter :: material_configFile = 'material.config' character(len=64), parameter :: material_configFile = 'material.config'
character(len=32), parameter :: material_partHomogenization = 'homogenization'
character(len=32), parameter :: material_partMicrostructure = 'microstructure'
character(len=32), parameter :: material_partPhase = 'phase'
character(len=32), parameter :: material_partTexture = 'texture'
!************************************* !*************************************
@ -68,13 +72,13 @@ subroutine material_init()
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
write(6,*) 'parsing homogenization...' write(6,*) 'parsing homogenization...'
call material_parseHomogenization(fileunit) call material_parseHomogenization(fileunit,material_partHomogenization)
write(6,*) 'parsing microstrcuture...' write(6,*) 'parsing microstrcuture...'
call material_parseMicrostructure(fileunit) call material_parseMicrostructure(fileunit,material_partMicrostructure)
write(6,*) 'parsing texture...' write(6,*) 'parsing texture...'
call material_parseTexture(fileunit) call material_parseTexture(fileunit,material_partTexture)
write(6,*) 'parsing phase...' write(6,*) 'parsing phase...'
call material_parsePhase(fileunit) call material_parsePhase(fileunit,material_partPhase)
close(fileunit) close(fileunit)
do i = 1,material_Nmicrostructure do i = 1,material_Nmicrostructure
@ -105,14 +109,14 @@ end subroutine
!********************************************************************* !*********************************************************************
subroutine material_parseHomogenization(file) subroutine material_parseHomogenization(file,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
character(len=32), parameter :: myPart = 'homogenization' character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 2 integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
@ -161,14 +165,14 @@ subroutine material_parseHomogenization(file)
!********************************************************************* !*********************************************************************
subroutine material_parseMicrostructure(file) subroutine material_parseMicrostructure(file,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
character(len=32), parameter :: myPart = 'microstructure' character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7 integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
@ -231,14 +235,14 @@ subroutine material_parseMicrostructure(file)
!********************************************************************* !*********************************************************************
subroutine material_parsePhase(file) subroutine material_parsePhase(file,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt use prec, only: pInt
use IO use IO
implicit none implicit none
character(len=32), parameter :: myPart = 'phase' character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 2 integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
@ -246,7 +250,7 @@ subroutine material_parsePhase(file)
character(len=64) tag character(len=64) tag
character(len=1024) line character(len=1024) line
Nsections= IO_countSections(file,myPart) Nsections = IO_countSections(file,myPart)
material_Nphase = Nsections material_Nphase = Nsections
allocate(phase_name(Nsections)); phase_name = '' allocate(phase_name(Nsections)); phase_name = ''
allocate(phase_constitution(Nsections)); phase_constitution = '' allocate(phase_constitution(Nsections)); phase_constitution = ''
@ -291,7 +295,7 @@ subroutine material_parsePhase(file)
!********************************************************************* !*********************************************************************
subroutine material_parseTexture(file) subroutine material_parseTexture(file,myPart)
!********************************************************************* !*********************************************************************
use prec, only: pInt, pReal use prec, only: pInt, pReal
@ -299,7 +303,7 @@ subroutine material_parseTexture(file)
use math, only: inRad use math, only: inRad
implicit none implicit none
character(len=32), parameter :: myPart = 'texture' character(len=*), intent(in) :: myPart
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 13 integer(pInt), parameter :: maxNchunks = 13
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions

View File

@ -1210,7 +1210,7 @@ end subroutine
real(pReal), intent(in) :: omega real(pReal), intent(in) :: omega
real(pReal), dimension(3) :: axisNrm real(pReal), dimension(3) :: axisNrm
real(pReal), dimension(3,3) :: math_RodrigToR real(pReal), dimension(3,3) :: math_RodrigToR
real(pReal) s,c real(pReal) s,c,c1
integer(pInt) i integer(pInt) i
forall (i=1:3) axisNrm(i) = axis(i)/dsqrt(math_mul3x3(axis,axis)) forall (i=1:3) axisNrm(i) = axis(i)/dsqrt(math_mul3x3(axis,axis))

View File

@ -33,8 +33,8 @@
include "IO.f90" ! uses prec, debug, math include "IO.f90" ! uses prec, debug, math
include "FEsolving.f90" ! uses prec, IO include "FEsolving.f90" ! uses prec, IO
include "mesh.f90" ! uses prec, IO, math, FEsolving include "mesh.f90" ! uses prec, IO, math, FEsolving
include "lattice.f90" ! uses prec, math
include "material.f90" ! uses prec, math, IO, mesh include "material.f90" ! uses prec, math, IO, mesh
include "lattice.f90" ! uses prec, math
include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug
@ -168,7 +168,9 @@
debug_InnerLoopDistribution = 0_pInt debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt debug_OuterLoopDistribution = 0_pInt
debug_cumLpTicks = 0 debug_cumLpTicks = 0
debug_cumDotStateTicks = 0
debug_cumLpCalls = 0_pInt debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt
endif endif
endif endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence if (cptim > theTime .or. theInc /= inc) then ! reached convergence

View File

@ -168,7 +168,9 @@
debug_InnerLoopDistribution = 0_pInt debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt debug_OuterLoopDistribution = 0_pInt
debug_cumLpTicks = 0 debug_cumLpTicks = 0
debug_cumDotStateTicks = 0
debug_cumLpCalls = 0_pInt debug_cumLpCalls = 0_pInt
debug_cumDotStateCalls = 0_pInt
endif endif
endif endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence if (cptim > theTime .or. theInc /= inc) then ! reached convergence