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

lattice: (re)introduced _symmetryType function to replace unsafe lookup array 

numerics: polishing

mpie_cpfem_marc: polishing

..powerlaw: aware of symmetryType function

crystallite: aware of symmetryType function, smaller leapfrog acceleration

IO: new warning 101

CPFEM: range of odd stress is now -1e15...+1e15, H_sym is used for stiffness
This commit is contained in:
Philip Eisenlohr 2010-11-03 14:58:11 +00:00
parent 763c20b302
commit 0dd99cb965
7 changed files with 97 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
code/CPFEM.f90
View File

@ -352,8 +352,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
if (selectiveDebugger) then
!$OMP CRITICAL (write2out)
write(6,'(a32,x,i8,x,i2,/,4(3(e20.8,x),/))') '°°° AGED state of grain 1, element ip',cp_en,IP, &
constitutive_state(1,IP,cp_en)%p
write(6,'(a,x,i8,x,i2,/,4(3(e20.8,x),/))') '<< cpfem >> AGED state of grain 1, element ip',&
cp_en,IP, constitutive_state(1,IP,cp_en)%p
!$OMP END CRITICAL (write2out)
endif
do k = 1,mesh_NcpElems
@ -433,16 +433,14 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
if (terminallyIll .or. outdatedFFN1 .or. any(abs(ffn1 - materialpoint_F(:,:,IP,cp_en)) > defgradTolerance)) then
if (.not. terminallyIll .and. .not. outdatedFFN1) then
!$OMP CRITICAL (write2out)
write(6,'(a32,x,i5,x,i2)') '°°° OUTDATED at element ip',cp_en,IP
write(6,'(a32,/,3(3(f10.6,x),/))') ' FFN was:',materialpoint_F(:,1,IP,cp_en), &
materialpoint_F(:,2,IP,cp_en), &
materialpoint_F(:,3,IP,cp_en)
write(6,'(a32,/,3(3(f10.6,x),/))') ' FFN1 now:',ffn1(:,1),ffn1(:,2),ffn1(:,3)
write(6,'(a,x,i5,x,i2)') '<< cpfem >> OUTDATED at element ip',cp_en,IP
write(6,'(a,/,3(3(f10.6,x),/))') ' FFN1 old:',materialpoint_F(1:3,:,IP,cp_en)
write(6,'(a,/,3(3(f10.6,x),/))') ' FFN1 now:',ffn1(1:3,:)
!$OMP END CRITICAL (write2out)
outdatedFFN1 = .true.
endif
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd
rnd = 2.0_pReal * rnd - 1.0_pReal
CPFEM_cs(:,IP,cp_en) = rnd*CPFEM_odd_stress
CPFEM_dcsde(:,:,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6)
@ -480,7 +478,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
if ( terminallyIll ) then
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd
rnd = 2.0_pReal * rnd - 1.0_pReal
CPFEM_cs(:,IP,cp_en) = rnd*CPFEM_odd_stress
CPFEM_dcsde(:,:,IP,cp_en) = CPFEM_odd_jacobian*math_identity2nd(6)
else
@ -492,7 +490,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! translate from dP/dF to dCS/dE
H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
! forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
H(i,j,k,l) = H(i,j,k,l) + &
materialpoint_F(j,m,IP,cp_en) * &
materialpoint_F(l,n,IP,cp_en) * &
@ -501,9 +498,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
0.5_pReal*(math_I3(i,k)*Kirchhoff(j,l) + math_I3(j,l)*Kirchhoff(i,k) + &
math_I3(i,l)*Kirchhoff(j,k) + math_I3(j,k)*Kirchhoff(i,l))
enddo; enddo; enddo; enddo; enddo; enddo
! forall(i=1:3,j=1:3,k=1:3,l=1:3) &
! H_sym(i,j,k,l) = 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k)) ! where to use this symmetric version ??
CPFEM_dcsde(:,:,IP,cp_en) = math_Mandel3333to66(J_inverse*H) ! should this use the symmetrized H ??
do i=1,3; do j=1,3; do k=1,3; do l=1,3
H_sym(i,j,k,l) = 0.25_pReal*(H(i,j,k,l)+H(j,i,k,l)+H(i,j,l,k)+H(j,i,l,k))
enddo; enddo; enddo; enddo
CPFEM_dcsde(:,:,IP,cp_en) = math_Mandel3333to66(J_inverse*H_sym)
endif
endif
@ -515,7 +513,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
CPFEM_dcsde = CPFEM_dcsde_knownGood ! --+>> RESTORE CONSISTENT JACOBIAN FROM FORMER CONVERGED INC
end if
call random_number(rnd)
if (rnd < 0.5_pReal) rnd = 1.0_pReal - rnd
rnd = 2.0_pReal * rnd - 1.0_pReal
materialpoint_Temperature(IP,cp_en) = Temperature
materialpoint_F0(:,:,IP,cp_en) = ffn
materialpoint_F(:,:,IP,cp_en) = ffn1
@ -539,8 +537,9 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! copy P and dPdF to the output variables
pstress(:,:) = materialpoint_P(:,:,IP,cp_en)
dPdF(:,:,:,:) = materialpoint_dPdF(:,:,:,:,IP,cp_en)
if ((debugger .and. selectiveDebugger) .and. &
mode < 6) then
selectiveDebugger = (cp_en == debug_e .and. IP == debug_i)
if (selectiveDebugger .and. mode < 6) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,i2,x,a,x,i4,/,6(f10.3,x)/)') 'stress/MPa at ip', IP, 'el', cp_en, cauchyStress/1e6
write(6,'(a,x,i2,x,a,x,i4,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip', IP, 'el', cp_en, jacobian/1e9

8
code/IO.f90
View File

@ -1373,14 +1373,16 @@ endfunction
character(len=80) msg
select case (ID)
case (101)
msg = '+ crystallite debugging off... +'
case (600)
msg = '+ Crystallite responds elastically +'
msg = '+ crystallite responds elastically +'
case (650)
msg = '+ Polar decomposition failed +'
msg = '+ polar decomposition failed +'
case (700)
msg = '+ unknown crystal symmetry +'
case default
msg = '+ Unknown warning number... +'
msg = '+ unknown warning number... +'
end select
!$OMP CRITICAL (write2out)

19
code/constitutive_phenopowerlaw.f90
View File

@ -29,12 +29,12 @@
!c44 46.7e9
!
!gdot0_slip 0.001
!1_by_m_slip 50
!n_slip 50
!tau0_slip 65e6 22e6 52e6 50e6 # per family
!tausat_slip 80e6 180e6 140e6 140e6 # per family
!w0_slip 1
!gdot0_twin 0.001
!1_by_m_twin 50
!n_twin 50
!tau0_twin 52e6 52e6 52e6 52e6 # per family
!s_pr 50e6 # push-up stress for slip saturation due to twinning
!twin_b 2
@ -134,7 +134,8 @@ subroutine constitutive_phenopowerlaw_init(file)
use IO
use material
use lattice, only: lattice_initializeStructure, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
use lattice, only: lattice_initializeStructure, lattice_symmetryType, &
lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_maxNinteraction, lattice_NslipSystem, lattice_NtwinSystem, &
lattice_interactionSlipSlip, &
lattice_interactionSlipTwin, &
@ -400,15 +401,17 @@ subroutine constitutive_phenopowerlaw_init(file)
constitutive_phenopowerlaw_sizeState(i) = constitutive_phenopowerlaw_totalNslip(i)+ &
constitutive_phenopowerlaw_totalNtwin(i)+ 2 ! s_slip, s_twin, sum(gamma), sum(f)
select case (constitutive_phenopowerlaw_structure(i)) ! assign elasticity tensor
case(1:2) ! cubic(s)
myStructure = constitutive_phenopowerlaw_structure(i)
select case (lattice_symmetryType(myStructure)) ! assign elasticity tensor
case(1) ! cubic(s)
forall(k=1:3)
forall(j=1:3) &
constitutive_phenopowerlaw_Cslip_66(k,j,i) = constitutive_phenopowerlaw_C12(i)
constitutive_phenopowerlaw_Cslip_66(k,k,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(k+3,k+3,i) = constitutive_phenopowerlaw_C44(i)
end forall
case(3:) ! hex
case(2) ! hex
constitutive_phenopowerlaw_Cslip_66(1,1,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(2,2,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(3,3,i) = constitutive_phenopowerlaw_C33(i)
@ -426,8 +429,6 @@ subroutine constitutive_phenopowerlaw_init(file)
constitutive_phenopowerlaw_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_phenopowerlaw_Cslip_66(:,:,i)))
myStructure = constitutive_phenopowerlaw_structure(i)
do f = 1,lattice_maxNslipFamily ! >>> interaction slip -- X
index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1,i))
do j = 1,constitutive_phenopowerlaw_Nslip(f,i) ! loop over (active) systems in my family (slip)
@ -479,6 +480,8 @@ subroutine constitutive_phenopowerlaw_init(file)
enddo; enddo
! report to out file...
enddo
return

44
code/crystallite.f90
View File

@ -93,7 +93,11 @@ subroutine crystallite_init(Temperature)
subStepSizeCryst, &
stepIncreaseCryst
use math, only: math_I3, &
math_EulerToR
math_EulerToR, &
math_inv3x3, &
math_transpose3x3, &
math_mul33xx33, &
math_mul33x33
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
use mesh, only: mesh_element, &
@ -102,9 +106,13 @@ subroutine crystallite_init(Temperature)
mesh_maxNipNeighbors
use IO
use material
use lattice, only: lattice_symmetryTypes
use lattice, only: lattice_symmetryType, &
lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem
use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_label, &
constitutive_phenopowerlaw_structure
constitutive_phenopowerlaw_structure, &
constitutive_phenopowerlaw_Nslip
use constitutive_titanmod, only: constitutive_titanmod_label, &
constitutive_titanmod_structure
use constitutive_dislotwin, only: constitutive_dislotwin_label, &
@ -123,6 +131,8 @@ subroutine crystallite_init(Temperature)
!*** local variables ***!
integer(pInt), parameter :: maxNchunks = 2
integer(pInt), dimension(1+2*maxNchunks) :: positions
character(len=64) tag
character(len=1024) line
integer(pInt) g, & ! grain number
i, & ! integration point number
e, & ! element number
@ -131,12 +141,11 @@ subroutine crystallite_init(Temperature)
eMax, & ! maximum number of elements
nMax, & ! maximum number of ip neighbors
myNgrains, & ! number of grains in current IP
myCrystallite ! crystallite of current elem
integer(pInt) section, j,p, output, mySize
character(len=64) tag
character(len=1024) line
integer(pInt) myStructure, & ! lattice structure
myPhase
myCrystallite, & ! crystallite of current elem
section, f,j,k,p, output, mySize, index_myFamily, &
myStructure, & ! lattice structure
myPhase, &
myMat
write(6,*)
write(6,*) '<<<+- crystallite init -+>>>'
@ -281,7 +290,7 @@ subroutine crystallite_init(Temperature)
do g = 1,myNgrains
crystallite_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption
crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
crystallite_Fe(:,:,g,i,e) = transpose(crystallite_Fp0(:,:,g,i,e))
crystallite_Fe(:,:,g,i,e) = math_transpose3x3(crystallite_Fp0(:,:,g,i,e))
crystallite_F0(:,:,g,i,e) = math_I3
crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e)
crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
@ -299,20 +308,21 @@ subroutine crystallite_init(Temperature)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e))
myPhase = material_phase(g,i,e)
myMat = phase_constitutionInstance(myPhase)
select case (phase_constitution(myPhase))
case (constitutive_phenopowerlaw_label)
myStructure = constitutive_phenopowerlaw_structure(phase_constitutionInstance(myPhase))
myStructure = constitutive_phenopowerlaw_structure(myMat)
case (constitutive_titanmod_label)
myStructure = constitutive_titanmod_structure(phase_constitutionInstance(myPhase))
myStructure = constitutive_titanmod_structure(myMat)
case (constitutive_dislotwin_label)
myStructure = constitutive_dislotwin_structure(phase_constitutionInstance(myPhase))
myStructure = constitutive_dislotwin_structure(myMat)
case (constitutive_nonlocal_label)
myStructure = constitutive_nonlocal_structure(phase_constitutionInstance(myPhase))
myStructure = constitutive_nonlocal_structure(myMat)
case default
myStructure = -1_pInt ! does this happen for j2 material?
end select
if (myStructure>0_pInt) then
crystallite_symmetryID(g,i,e)=lattice_symmetryTypes(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3(hex)=>2
crystallite_symmetryID(g,i,e) = lattice_symmetryType(myStructure) ! structure = 1(fcc) or 2(bcc) => 1; 3(hex)=>2
endif
enddo
enddo
@ -2386,7 +2396,7 @@ endfunction
! start LpLoop with no acceleration
NiterationStress = 0_pInt
leapfrog = 1.0_pReal
maxleap = 1024.0_pReal
maxleap = 16.0_pReal
jacoCounter = 0_pInt
LpLoop: do
@ -2523,7 +2533,7 @@ LpLoop: do
Lpguess_old = Lpguess
! accelerate?
if (NiterationStress > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog
if (NiterationStress > 1 .and. leapfrog < maxleap) leapfrog = leapfrog + 1.0_pReal
endif
! leapfrog to updated Lp

29
code/lattice.f90
View File

@ -26,9 +26,6 @@ integer(pInt), parameter :: lattice_maxNslip = 48 ! max # of sli
integer(pInt), parameter :: lattice_maxNtwin = 24 ! max # of twin systems over lattice structures
integer(pInt), parameter :: lattice_maxNinteraction = 20 ! max # of interaction types (in hardening matrix part)
integer(pInt), parameter, dimension(3) :: lattice_symmetryTypes =(/1, 1, 2/) ! maps crystal structures to symmetry tpyes
integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, &
interactionSlipTwin, &
interactionTwinSlip, &
@ -661,6 +658,30 @@ CONTAINS
!* - lattice_initializeStructure
!****************************************
pure function lattice_symmetryType(structID)
!**************************************
!* maps structure to symmetry type *
!* fcc(1) and bcc(2) are cubic(1) *
!* hex(3+) is hexagonal(2) *
!**************************************
implicit none
integer(pInt), intent(in) :: structID
integer(pInt) lattice_symmetryType
select case(structID)
case (1,2)
lattice_symmetryType = 1_pInt
case (3:)
lattice_symmetryType = 2_pInt
case default
lattice_symmetryType = 0_pInt
end select
return
end function
subroutine lattice_init()
!**************************************
@ -684,7 +705,7 @@ subroutine lattice_init()
! lattice_Nstructure = Nsections + 2_pInt ! most conservative assumption
close(fileunit)
write(6,'(a16,x,i5)') '# sections:',Nsections
write(6,'(a16,x,i5)') '# phases:',Nsections
write(6,'(a16,x,i5)') '# structures:',lattice_Nstructure
write(6,*)

6
code/mpie_cpfem_marc.f90
View File

@ -97,14 +97,14 @@ END MODULE
include "mesh.f90" ! uses prec, math, IO, FEsolving
include "material.f90" ! uses prec, math, IO, mesh
include "lattice.f90" ! uses prec, math, IO, material
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_titanmod.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_dislotwin.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_nonlocal.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug
include "crystallite.f90" ! uses prec, math, IO, numerics
include "homogenization_isostrain.f90" ! uses prec, math, IO,
include "crystallite.f90" ! uses prec, math, IO, numerics, Fesolving, material, mesh, constitutive
include "homogenization_isostrain.f90" ! uses prec, math, IO
include "homogenization_RGC.f90" ! uses prec, math, IO, numerics, mesh: added <<<updated 31.07.2009>>>
include "homogenization.f90" ! uses prec, math, IO, numerics
include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite

4
code/numerics.f90
View File

@ -110,8 +110,8 @@ subroutine numerics_init()
rTol_crystalliteTemperature = 1.0e-6_pReal
rTol_crystalliteStress = 1.0e-6_pReal
aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here)
integrator = 1
integratorStiffness = 1
integrator = 1 ! fix-point iteration
integratorStiffness = 1 ! fix-point iteration
!* RGC parameters: added <<<updated 17.12.2009>>> with moderate setting
absTol_RGC = 1.0e+4