diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index a80adcacf..74033fa46 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -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 diff --git a/code/IO.f90 b/code/IO.f90 index 2b46db842..3202f1772 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -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) diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 79020d14e..a00f21757 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -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 @@ -592,7 +595,7 @@ subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - + endsubroutine diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 5c5afb197..a5b624545 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -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,15 +106,19 @@ 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 - use constitutive_titanmod, only: constitutive_titanmod_label, & - constitutive_titanmod_structure + constitutive_phenopowerlaw_structure, & + constitutive_phenopowerlaw_Nslip + use constitutive_titanmod, only: constitutive_titanmod_label, & + constitutive_titanmod_structure use constitutive_dislotwin, only: constitutive_dislotwin_label, & - constitutive_dislotwin_structure + constitutive_dislotwin_structure use constitutive_nonlocal, only: constitutive_nonlocal_label, & - constitutive_nonlocal_structure + constitutive_nonlocal_structure implicit none integer(pInt), parameter :: file = 200 @@ -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 diff --git a/code/lattice.f90 b/code/lattice.f90 index e446ac914..6061906f6 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -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,*) diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index 912e5df22..a0edc9057 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -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_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_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug + include "constitutive_titanmod.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 <<>> include "homogenization.f90" ! uses prec, math, IO, numerics include "CPFEM.f90" ! uses prec, math, IO, numerics, debug, FEsolving, mesh, lattice, constitutive, crystallite diff --git a/code/numerics.f90 b/code/numerics.f90 index 2c7f1c706..9647934ef 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -15,8 +15,8 @@ integer(pInt) iJacoStiffness, & ! freque nState, & ! state loop limit nStress, & ! stress loop limit pert_method, & ! method used in perturbation technique for tangent - integrator, & ! method used for state integration - integratorStiffness ! method used for stiffness state integration + integrator, & ! method used for state integration + integratorStiffness ! method used for stiffness state integration real(pReal) relevantStrain, & ! strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) defgradTolerance, & ! deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) pert_Fg, & ! strain perturbation for FEM Jacobi @@ -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 <<>> with moderate setting absTol_RGC = 1.0e+4