From 5aff0d96b6369408c5ebed9721f5699aa42ac0f6 Mon Sep 17 00:00:00 2001 From: achalhp Date: Wed, 18 Oct 2023 18:07:54 +0530 Subject: [PATCH] 18-Oct-23-1 --- ...phase_mechanical_plastic_phenopowerlaw.f90 | 446 +++++++++++++++--- 1 file changed, 383 insertions(+), 63 deletions(-) diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index f0dc04869..df70b55f7 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -1,3 +1,17 @@ +! Copyright 2011-2022 Max-Planck-Institut für Eisenforschung GmbH +! +! DAMASK is free software: you can redistribute it and/or modify +! it under the terms of the GNU Affero General Public License as +! published by the Free Software Foundation, either version 3 of the +! License, or (at your option) any later version. +! +! This program is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Affero General Public License for more details. +! +! You should have received a copy of the GNU Affero General Public License +! along with this program. If not, see . !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH @@ -7,132 +21,158 @@ submodule(phase:plastic) phenopowerlaw type :: tParameters - real(pReal) :: & + real(pReal) :: & !< declaring real number variables dot_gamma_0_sl = 1.0_pReal, & !< reference shear strain rate for slip dot_gamma_0_tw = 1.0_pReal, & !< reference shear strain rate for twin n_sl = 1.0_pReal, & !< stress exponent for slip n_tw = 1.0_pReal, & !< stress exponent for twin f_sat_sl_tw = 1.0_pReal, & !< push-up factor for slip saturation due to twinning - c_1 = 1.0_pReal, & + c_1 = 1.0_pReal, & !< hardening parameters c_2 = 1.0_pReal, & c_3 = 1.0_pReal, & c_4 = 1.0_pReal, & h_0_sl_sl = 1.0_pReal, & !< reference hardening slip - slip h_0_tw_sl = 1.0_pReal, & !< reference hardening twin - slip h_0_tw_tw = 1.0_pReal, & !< reference hardening twin - twin - a_sl = 1.0_pReal - real(pReal), allocatable, dimension(:) :: & + h_0_tw_tw_nuc = 1.0_pReal, & !< Achal twin nucleation + h_0_tw_tw_grt = 1.0_pReal, & !< Achal twin growth + a_sl = 1.0_pReal, & !< non-Schmid Coefficient + chkstep_nucl = 1.0_pReal, & !< Achal Monte Carlo sampling frequency (for twin nucleation) + chkstep_grow = 1.0_pReal, & !< Achal Monte Carlo sampling frequency (for twin growth) + chkgrowth_twin = 1.0_pReal, & !< Achal flag for twin growth + prefdecay_slip = 1.0_pReal, & !< Achal pre-factor for the flag_slip decay + twin_inclusion = 1.0_pReal !< Achal flag_slip to introduce specific twin vol fraction to particular phase + + real(pReal), allocatable, dimension(:) :: & !< 1D array xi_inf_sl, & !< maximum critical shear stress for slip h_int, & !< per family hardening activity (optional) gamma_char !< characteristic shear for twins - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & !< 2D array h_sl_sl, & !< slip resistance from slip activity h_sl_tw, & !< slip resistance from twin activity h_tw_sl, & !< twin resistance from slip activity - h_tw_tw !< twin resistance from twin activity - real(pReal), allocatable, dimension(:,:,:) :: & - P_sl, & - P_tw, & - P_nS_pos, & - P_nS_neg + !h_tw_tw!, & !< twin resistance from twin activity + h_tw_tw_grow, & + h_tw_tw_nucl + real(pReal), allocatable, dimension(:,:,:) :: & !< 3D array + P_sl, & !< Schmid slip + P_tw, & !< Schmid twin + P_nS_pos, & !< Non Schmid +ve + P_nS_neg, & !< Non Schmid -ve + Corrs_Matr !< Achal Correspondence Matrix integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw !< total number of active twin systems - logical :: & + logical :: & !< flag if non-Schmid slip is Active nonSchmidActive = .false. - character(len=pStringLen), allocatable, dimension(:) :: & + character(len=pStringLen), allocatable, dimension(:) :: & !< Storing output messages output - character(len=:), allocatable, dimension(:) :: & + character(len=:), allocatable, dimension(:) :: & !< Labels of active slip and twin systems systems_sl, & systems_tw end type tParameters - type :: tIndexDotState + type :: tIndexDotState !< Index to access variables in 2D tPhenopowerlawState integer, dimension(2) :: & xi_sl, & - xi_tw, & + !xi_tw, & gamma_sl, & - gamma_tw + gamma_tw, & + xi_tw_nucl, & + xi_tw_grow end type tIndexDotState - type :: tPhenopowerlawState + type :: tPhenopowerlawState !< state variables at material points real(pReal), pointer, dimension(:,:) :: & - xi_sl, & - xi_tw, & - gamma_sl, & - gamma_tw + xi_sl, & !< critical shear stress for slip + !xi_tw, & !< critical shear stress for twin + gamma_sl, & !< shear strain due to slip + gamma_tw, & !< shear strain due to twin + xi_tw_nucl, & + xi_tw_grow, & + f_tw_nucl, & + f_tw_grow, & + fmc_tw_nucl, & + fmc_tw_grwo + !real(pReal), pointer, dimension(:) :: & + ! variant_twin, & !< flag used to assign twin variant + ! frozen !< 0 to 1 + end type tPhenopowerlawState !-------------------------------------------------------------------------------------------------- ! containers for parameters, dot state index, and state +! containers are arrays used to store parameters, state indices and state variables for multiple material points type(tParameters), allocatable, dimension(:) :: param type(tIndexDotState), allocatable, dimension(:) :: indexDotState type(tPhenopowerlawState), allocatable, dimension(:) :: state contains +!< "contains" means the above types are used in below functions. !-------------------------------------------------------------------------------------------------- !> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function plastic_phenopowerlaw_init() result(myPlasticity) +module function plastic_phenopowerlaw_init() result(myPlasticity) !< function "plastic_phenopowerlaw_init" returns myPlasticity - logical, dimension(:), allocatable :: myPlasticity - integer :: & + logical, dimension(:), allocatable :: myPlasticity !< logical array: phenopowerlaw is active or not + integer :: & !< integer variables to iterate ph, i, & Nmembers, & sizeState, sizeDotState, & startIndex, endIndex - integer, dimension(:), allocatable :: & + integer, dimension(:), allocatable :: & !< active number of slip and twin systems N_sl, N_tw real(pReal), dimension(:), allocatable :: & xi_0_sl, & !< initial critical shear stress for slip xi_0_tw, & !< initial critical shear stress for twin a !< non-Schmid coefficients - character(len=pStringLen) :: & + character(len=pStringLen) :: & !< to store error messages extmsg = '' - class(tNode), pointer :: & + class(tNode), pointer :: & !< pointers variables, tNode objects phases, & phase, & - mech, & - pl + mech, & !< mechanical model info + pl !< plastic model info - myPlasticity = plastic_active('phenopowerlaw') - if(count(myPlasticity) == 0) return + myPlasticity = plastic_active('phenopowerlaw') !< function defined in phase_mechanical_plastic + if(count(myPlasticity) == 0) return !< if phenopowerlaw is active or not print'(/,1x,a)', '<<<+- phase:mechanical:plastic:phenopowerlaw init -+>>>' - print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) + print'(/,a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT) !< print for which phases phenopowerlaw is active - phases => config_material%get('phase') + phases => config_material%get('phase') !< get "phase" objects from config_material object allocate(param(phases%length)) allocate(indexDotState(phases%length)) allocate(state(phases%length)) + - do ph = 1, phases%length + do ph = 1, phases%length !< looping for each phase if (.not. myPlasticity(ph)) cycle - associate(prm => param(ph), stt => state(ph), & + associate(prm => param(ph), stt => state(ph), & !< parameters, state var, dot state indices of current phase idx_dot => indexDotState(ph)) - phase => phases%get(ph) + phase => phases%get(ph) !< getting objects mech => phase%get('mechanical') pl => mech%get('plastic') !-------------------------------------------------------------------------------------------------- ! slip related parameters - N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) - prm%sum_N_sl = sum(abs(N_sl)) - slipActive: if (prm%sum_N_sl > 0) then - prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) + N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) !< Read and Store no. of slip systems + prm%sum_N_sl = sum(abs(N_sl)) !< total no. of slip systems + slipActive: if (prm%sum_N_sl > 0) then !< if +ve, active slip system exists + prm%systems_sl = lattice_labels_slip(N_sl,phase_lattice(ph)) !< prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph)) - if (phase_lattice(ph) == 'cI') then - a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray) + if (phase_lattice(ph) == 'cI') then !< cI: BCC + a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray) !< a: non-Schmid coefficients for projections if(size(a) > 0) prm%nonSchmidActive = .true. - prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) + prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1) !< function returns non-schmid projections prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1) else prm%P_nS_pos = prm%P_sl @@ -177,6 +217,14 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) prm%P_tw = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%h_tw_tw = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'),phase_lattice(ph)) prm%gamma_char = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),phase_cOverA(ph)) + !prm%h_tw_tw_nucl = + !prm%h_tw_tw_grow = + !prm%chkstep_nucl = + !prm%chkstep_grow = + !prm%chkgrowth_twin = + !prm%twin_inclusion = + !prm%prefdecay_slip = + xi_0_tw = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw)) @@ -200,6 +248,8 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) xi_0_tw = emptyRealArray allocate(prm%gamma_char,source=emptyRealArray) allocate(prm%h_tw_tw(0,0)) + !allocate(prm%h_tw_tw_nucl(0,0)) + !allocate(prm%h_tw_tw_grow(0,0)) end if twinActive !-------------------------------------------------------------------------------------------------- @@ -229,8 +279,16 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) ! allocate state arrays Nmembers = count(material_phaseID == ph) sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & - + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw - sizeState = sizeDotState + + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw !& + + size(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw & ! Why not size(['xi_tw_nucl','gamma_tw'])? + !+ size(['f_tw_nucl','f_tw_grow']) * prm%sum_N_tw & + !+ size(['variant_twin','frozen']) * prm%sum_N_tw & + sizeState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw !& + + size(['xi_tw_nucl','xi_tw_grow']) * prm%sum_N_tw & ! Why not size(['xi_tw_nucl','gamma_tw'])? + !+ size(['f_tw_nucl','f_tw_grow']) * prm%sum_N_tw & + !+ size(['fmc_tw_nucl','fmc_tw_grow']) * prm%sum_N_tw & + !+ size(['variant_twin','frozen']) * prm%sum_N_tw & call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely @@ -245,12 +303,12 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' - startIndex = endIndex + 1 - endIndex = endIndex + prm%sum_N_tw - idx_dot%xi_tw = [startIndex,endIndex] - stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:) - stt%xi_tw = spread(xi_0_tw, 2, Nmembers) - plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) + !startIndex = endIndex + 1 + !endIndex = endIndex + prm%sum_N_tw + !idx_dot%xi_tw = [startIndex,endIndex] + !stt%xi_tw => plasticState(ph)%state(startIndex:endIndex,:) + !stt%xi_tw = spread(xi_0_tw, 2, Nmembers) + !plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl @@ -265,6 +323,50 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) stt%gamma_tw => plasticState(ph)%state(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_tw + idx_dot%xi_tw_nucl = [startIndex,endIndex] + stt%xi_tw_nucl => plasticState(ph)%state(startIndex:endIndex,:) + stt%xi_tw_nucl = spread(xi_0_tw, 2, Nmembers) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%sum_N_tw + idx_dot%xi_tw_grow = [startIndex,endIndex] + stt%xi_tw_grow => plasticState(ph)%state(startIndex:endIndex,:) + stt%xi_tw_grow = spread(xi_0_tw, 2, Nmembers) + plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + + !startIndex = endIndex + 1 + !endIndex = endIndex + prm%sum_N_tw + !idx_dot%f_twin_nucl = [startIndex,endIndex] + !stt%f_twin_nucl => plasticState(ph)%state(startIndex:endIndex,:) + !plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + + !startIndex = endIndex + 1 + !endIndex = endIndex + prm%sum_N_tw + !idx_dot%f_twin_grow = [startIndex,endIndex] + !stt%f_twin_grow => plasticState(ph)%state(startIndex:endIndex,:) + !plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + + !startIndex = endIndex + 1 + !endIndex = endIndex + prm%sum_N_tw + !idx_dot%fmc_twin_nucl = [startIndex,endIndex] + !stt%fmc_twin_nucl => plasticState(ph)%state(startIndex:endIndex,:) + !plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + + !startIndex = endIndex + 1 + !endIndex = endIndex + prm%sum_N_tw + !idx_dot%fmc_twin_grow = [startIndex,endIndex] + !stt%fmc_twin_grow => plasticState(ph)%state(startIndex:endIndex,:) + !plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + + !startIndex = endIndex + 1 + !endIndex = endIndex + prm%sum_N_tw + !idx_dot%frozen = [startIndex,endIndex] + !stt%frozen => plasticState(ph)%state(startIndex:endIndex,:) + !plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) + end associate !-------------------------------------------------------------------------------------------------- @@ -299,9 +401,9 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_sl_pos,dot_gamma_sl_neg, & ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg + real(pReal), dimension(param(ph)%sum_N_tw) :: & - dot_gamma_tw,ddot_gamma_dtau_tw - + dot_gamma_tw,ddot_gamma_dtau_tw Lp = 0.0_pReal dLp_dMp = 0.0_pReal @@ -332,7 +434,7 @@ end subroutine phenopowerlaw_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- -module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) +module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) !< Rate of change of state variables real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -348,13 +450,23 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) real(pReal), dimension(param(ph)%sum_N_sl) :: & dot_gamma_sl_pos,dot_gamma_sl_neg, & right_SlipSlip + !real(pReal), dimension(param(ph)%sum_N_tw) :: & + !fdot_twin_nucl, fdot_twin_grow associate(prm => param(ph), stt => state(ph), & dot_xi_sl => dotState(indexDotState(ph)%xi_sl(1):indexDotState(ph)%xi_sl(2)), & - dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), & + !dot_xi_tw => dotState(indexDotState(ph)%xi_tw(1):indexDotState(ph)%xi_tw(2)), & dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), & - dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2))) + !dot_gamma_tw => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2))!, & + dot_gamma_tw_nucl => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)), & + dot_gamma_tw_grow => dotState(indexDotState(ph)%gamma_tw(1):indexDotState(ph)%gamma_tw(2)), & + dot_xi_tw_nucl => dotState(indexDotState(ph)%xi_tw_nucl(1):indexDotState(ph)%xi_tw_nucl(2)), & + dot_xi_tw_grow => dotState(indexDotState(ph)%xi_tw_grow(1):indexDotState(ph)%xi_tw_grow(2)) + + !sumGamma + !sumF_nucl + !sumF_grow call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg) dot_gamma_sl = abs(dot_gamma_sl_pos+dot_gamma_sl_neg) @@ -368,15 +480,219 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) dot_xi_sl = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & * matmul(prm%h_sl_sl,dot_gamma_sl*right_SlipSlip) & + matmul(prm%h_sl_tw,dot_gamma_tw) +!< Rate of change of critical shear stress + !dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & + ! * matmul(prm%h_tw_sl,dot_gamma_sl) & + ! + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw) - dot_xi_tw = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & - * matmul(prm%h_tw_sl,dot_gamma_sl) & - + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot_gamma_tw) + dot_xi_tw_nucl = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & + * matmul(prm%h_tw_sl,dot_gamma_sl) & + + prm%h_0_tw_tw_nucl * sumF_nucl**prm%c_4 * matmul(prm%h_tw_tw_nucl,dot_gamma_tw) + + dot_xi_tw_grow = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & + * matmul(prm%h_tw_sl,dot_gamma_sl) & + + prm%h_0_tw_tw_grow * sumF_grow**prm%c_4 * matmul(prm%h_tw_tw_grow,dot_gamma_tw) end associate end function phenopowerlaw_dotState +!-------------------------------------------------------------------------------------------------- +!> @brief calculates instantaneous incremental change of kinematics and associated jump state +!-------------------------------------------------------------------------------------------------- +! module subroutine plastic_kinematic_deltaFp(twin) +! use prec, only: & +! dNeq, & +! dEq0 +! #ifdef DEBUG +! use debug, only: & +! debug_level, & +! debug_constitutive,& +! debug_levelExtensive, & +! debug_levelSelective +! #endif + +! use mesh, only: & +! mesh_element, & +! mesh_ipNeighborhood, & +! mesh_ipCoordinates, & +! mesh_ipVolume, & +! mesh_ipAreaNormal, & +! mesh_ipArea, & +! FE_NipNeighbors, & +! mesh_maxNipNeighbors, & +! FE_geomtype, & +! FE_celltype + +! use lattice +! use math, only: & +! math_I3 + +! use material, only: & +! phaseAt, phasememberAt, & +! phase_plasticityInstance + +! implicit none +! integer(pInt) :: & +! ph, of, instance, & +! neighbor_el, & !< element number of neighboring material point +! neighbor_ip, & !< integration point of neighboring material point +! np, & !< neighbor phase +! no, n !< nieghbor offset and index for loop at neighbor + +! integer(pInt), intent(in) :: & +! el, & !< element index +! ip, & !< integration point index +! ipc !< grain index +! real(pReal), dimension(3,3), intent(out) :: & +! deltaFp +! logical , intent(out) :: & +! twinJump +! ! real(pReal), dimension(3,3,param(instance)%totalNslip) :: & +! ! CorrespondanceMatrix +! integer(pInt), dimension(52) :: & +! twin_el_incl +! real(pReal), dimension(6) :: & +! neighbor_stt +! real(pReal) :: & +! random, random1 +! integer(pInt) :: & +! i,j,var_growth,var_nucl +! var_growth = 0_pInt +! var_nucl = 0_pInt +! ph = phaseAt(ipc, ip, el) +! of = phasememberAt(ipc, ip, el) +! instance = phase_plasticityInstance(ph) + +! associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) + +! twinJump = .false. +! deltaFp = math_I3 + +! ! for eshelby circular inclusion +! twin_el_incl = (/ 10913,10914,10915,10916,10917,10918,10919,10920,10921,10922,10923,10924,10925, & +! 10993,10994,10995,10996,10997,10998,10999,11000,11001,11002,11074,11075,11076, & +! 11077,11078,11079,10751,10752,10753,10754,10755,10756,10757,10758,10759,10760, & +! 10670,10671,10672,10673,10674,10675,10676,10677,10678,10679,10680,10681,10682 /) +! ! TwinLooptest: do i=1_pInt, prm%totalNtwin +! ! write(6,*)'CorrespondenceMatrix for system',i, prm%CorrespondanceMatrix(:,:,i) +! ! enddo TwinLooptest + + +! !Saving the neighbor information in an array +! NeighborLoop1: do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! only 4 neighbors for quasi 2D (1 element in z direction) +! neighbor_el = mesh_ipNeighborhood(1,n,ip,el) +! neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) +! np = phaseAt(1,neighbor_ip,neighbor_el) +! no = phasememberAt(1,neighbor_ip,neighbor_el) +! neighbor_stt(n) = state(phase_plasticityInstance(np))%variant_twin(no) +! enddo NeighborLoop1 + +! !checking if any of my neighbor is twinned if yes recognize the variant and exit +! ! NeighborLoop2: do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! only 4 neighbors for quasi 2D (1 element in z direction) +! ! neighbor_el = mesh_ipNeighborhood(1,n,ip,el) +! ! neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) +! ! np = phaseAt(1,neighbor_ip,neighbor_el) +! ! ! if(of == 1) write(6,*)'phaseAt neighbor_ip of neighbor_el', np +! ! no = phasememberAt(1,neighbor_ip,neighbor_el) +! ! ! if(of == 1) write(6,*)'phasememberAt at neighbor_ip of neighbor_el', no +! ! if (state(phase_plasticityInstance(np))%variant_twin(no) > 0_pInt) then +! ! var_growth = state(phase_plasticityInstance(np))%variant_twin(no) +! ! +! ! exit NeighborLoop2 +! ! endif +! ! enddo NeighborLoop2 + +! call RANDOM_NUMBER(random) +! call RANDOM_NUMBER(random1) +! ! Sampling: if (var_growth > 0_pInt) then +! ! ! write(6,*)'I am sampling for growth with variant',var_growth +! ! Ability_Growth: if (stt%f_twin_grow(var_growth,of) > stt%fmc_twin_grow(var_growth,of) & +! ! + prm%checkstep_grow) then +! ! stt%fmc_twin_grow(var_growth,of) = stt%fmc_twin_grow(var_growth,of) & +! ! + prm%checkstep_grow +! ! Success_Growth: if (random <= stt%f_twin_grow(var_growth,of) .or. ALL(neighbor_stt > 0_pReal)) then +! ! write(6,*)'growth sampling is successful for elem',el +! ! twinJump = .true. +! ! deltaFp = prm%CorrespondanceMatrix(:,:,var_growth) +! ! dlt%f_twin_grow(:,of) = 0.0_pReal - stt%f_twin_grow(:,of) +! ! dlt%f_twin_nucl(:,of) = 0.0_pReal - stt%f_twin_nucl(:,of) +! ! dlt%fmc_twin_grow(:,of) = 0.0_pReal - stt%fmc_twin_grow(:,of) +! ! dlt%fmc_twin_nucl(:,of) = 0.0_pReal - stt%fmc_twin_nucl(:,of) +! ! dlt%frozen(of) = 1.0_pReal - stt%frozen(of) +! ! dlt%variant_twin(of) = var_growth - stt%variant_twin(of) +! ! endif Success_Growth +! ! endif Ability_Growth + +! ! elseif (var_growth == 0_pInt .and. prm%checkgrowth_twin > 0_pReal ) then +! if (var_growth == 0_pInt .and. prm%checkgrowth_twin > 0_pReal ) then +! var_nucl = maxloc(stt%f_twin_nucl(:,of), dim=1) +! ! write(6,*)'I am sampling for nucleation with variant',var_nucl,stt%f_twin_nucl(var_nucl,of) +! Ability_Nucleation: if (stt%f_twin_nucl(var_nucl,of) > stt%fmc_twin_nucl(var_nucl,of) & +! + prm%checkstep_nucl) then +! stt%fmc_twin_nucl(var_nucl,of) = stt%fmc_twin_nucl(var_nucl,of) & +! + prm%checkstep_nucl +! Success_Nucleation: if (random <= stt%f_twin_nucl(var_nucl,of) & +! .and. random1 <= 0.20) then +! write(6,*)'nucleation sampling is successful for elem',el +! twinJump = .true. +! deltaFp = prm%CorrespondanceMatrix(:,:,var_nucl) +! dlt%f_twin_nucl(:,of) = 0.0_pReal - stt%f_twin_nucl(:,of) +! dlt%f_twin_grow(:,of) = 0.0_pReal - stt%f_twin_grow(:,of) +! dlt%fmc_twin_nucl(:,of) = 0.0_pReal - stt%fmc_twin_nucl(:,of) +! dlt%fmc_twin_grow(:,of) = 0.0_pReal - stt%fmc_twin_grow(:,of) +! dlt%frozen(of) = 1.0_pReal - stt%frozen(of) +! dlt%variant_twin(of) = var_nucl - stt%variant_twin(of) +! endif Success_Nucleation +! endif Ability_Nucleation +! endif +! ! endif Sampling +! end associate + +! end subroutine plastic_kinematic_deltaFp + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates (instantaneous) incremental change of microstructure +!-------------------------------------------------------------------------------------------------- +!subroutine plastic_phenopowerlaw_deltaState(instance,of) +! use prec, only: & +! dNeq, & +! dEq0 +! #ifdef DEBUG +! use debug, only: & +! debug_level, & +! debug_constitutive,& +! debug_levelExtensive, & +! debug_levelSelective +! #endif + +! implicit none +! integer(pInt), intent(in) :: & +! instance, & +! of +! +! associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) + +! #ifdef DEBUG +! if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt & +! .and. (of == prm%of_debug & +! .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then +! write(6,'(a)') '======= phenopowerlaw delta state =======' +! ! write(6,*) sense,state(instance)%sense(:,of) +! endif +! #endif + +! !-------------------------------------------------------------------------------------------------- +! dlt%f_twin_nucl(:,of) = 0.0_pReal +! dlt%f_twin_grow(:,of) = 0.0_pReal +! dlt%fmc_twin_nucl(:,of) = 0.0_pReal +! dlt%fmc_twin_grow(:,of) = 0.0_pReal +! dlt%frozen(of) = 0.0_pReal +! dlt%variant_twin(of) = 0.0_pInt +! +! end associate + +!end subroutine plastic_phenopowerlaw_deltaState !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file. @@ -494,9 +810,10 @@ end subroutine kinetics_sl !> @details Derivatives are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end. +! Mp: Mandel Stress, ph: , en: , dot_gamma_tw: , ddot_gamma_dtau_tw: !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_tw(Mp,ph,en,& - dot_gamma_tw,ddot_gamma_dtau_tw) + dot_gamma_tw,ddot_gamma_dtau_tw) !fdot_twin_nucl, fdot_twin_grow) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -505,7 +822,7 @@ pure subroutine kinetics_tw(Mp,ph,en,& en real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & - dot_gamma_tw + dot_gamma_tw !fdot_twin_nucl, fdot_twin_grow real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: & ddot_gamma_dtau_tw @@ -518,9 +835,11 @@ pure subroutine kinetics_tw(Mp,ph,en,& tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)] - where(tau_tw > 0.0_pReal) + where(tau_tw > 0.0_pReal) !and stt%frozen(en) < 0.9_pReal) dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw + !fdot_twin_nucl = + !fdot_twin_nucl = else where dot_gamma_tw = 0.0_pReal end where @@ -538,3 +857,4 @@ pure subroutine kinetics_tw(Mp,ph,en,& end subroutine kinetics_tw end submodule phenopowerlaw +