From f337847f351589433eae40eff808ce0537afc6e9 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 22 Jul 2009 16:07:19 +0000 Subject: [PATCH] quite some changes: # non-greedy memory allocation # generation of outputConstitutive to allow for script-based T16 extraction # exchange of phenomenological by more general phenopowerlaw # lattice is based on slip and twin families which can be treated as individual entities (switched on/off, separate hardening, etc.) # nicer debugging output # changed some error/warning codes # plus potentially some minor additional brushes here and there --- code/CPFEM.f90 | 46 +- code/IO.f90 | 60 +- code/constitutive.f90 | 111 ++- code/constitutive_dislobased.f90 | 817 +++++++++------------- code/constitutive_j2.f90 | 164 ++--- code/constitutive_phenomenological.f90 | 519 -------------- code/constitutive_phenopowerlaw.f90 | 912 +++++++++++++++++++++++++ code/crystallite.f90 | 473 +++++++------ code/debug.f90 | 14 +- code/homogenization.f90 | 151 ++-- code/lattice.f90 | 384 ++++++----- code/material.config | 106 +-- code/material.f90 | 40 +- code/mpie_cpfem_marc.f90 | 2 +- code/numerics.f90 | 33 +- 15 files changed, 2050 insertions(+), 1782 deletions(-) delete mode 100644 code/constitutive_phenomenological.f90 create mode 100644 code/constitutive_phenopowerlaw.f90 diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 830a20388..0386a7e0c 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -36,16 +36,16 @@ subroutine CPFEM_init() ! initialize stress and jacobian to zero allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal - allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde = 0.0_pReal - allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal + allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE = 0.0_pReal + allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE_knownGood = 0.0_pReal !$OMP CRITICAL (write2out) write(6,*) write(6,*) '<<<+- cpfem init -+>>>' write(6,*) write(6,'(a32,x,6(i5,x))') 'CPFEM_cs: ', shape(CPFEM_cs) - write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde: ', shape(CPFEM_dcsde) - write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood) + write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE) + write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood) write(6,*) write(6,*) 'parallelExecution: ', parallelExecution write(6,*) 'symmetricSolver: ', symmetricSolver @@ -61,7 +61,7 @@ endsubroutine !*** call the actual material model *** !*********************************************************************** subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian, ngens) - ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/de + ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE !*** variables and functions from other modules ***! use prec, only: pReal, & @@ -90,12 +90,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt use mesh, only: mesh_init, & mesh_FEasCP, & mesh_NcpElems, & - mesh_maxNips, & - mesh_element + mesh_maxNips use lattice, only: lattice_init use material, only: material_init, & - homogenization_maxNgrains, & - homogenization_Ngrains + homogenization_maxNgrains use constitutive, only: constitutive_init,& constitutive_state0,constitutive_state use crystallite, only: crystallite_init, & @@ -144,8 +142,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt real(pReal), dimension (3,3,3,3) :: H, & H_sym integer(pInt) cp_en, & ! crystal plasticity element number - e, & - g, & i, & j, & k, & @@ -183,6 +179,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt cp_en = mesh_FEasCP('elem',element) if (cp_en == 1 .and. IP == 1) then + write(6,*) write(6,*) '#####################################' write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') & 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& @@ -199,19 +196,16 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt if (mode == 1) then crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation - crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity - do e = 1,mesh_NcpElems - do g = 1,homogenization_Ngrains(mesh_element(3,e)) - do i = 1,mesh_maxNips - constitutive_state0(g,i,e)%p = constitutive_state(g,i,e)%p ! microstructure of crystallites - enddo - enddo - enddo - write(6,'(a10,/,4(3(f10.3,x),/))') 'aged state',constitutive_state(1,1,1)%p/1e6 - do e = 1,mesh_NcpElems - do i = 1,mesh_maxNips - if (homogenization_sizeState(i,e) > 0_pInt) & - homogenization_state0(i,e)%p = homogenization_state(i,e)%p ! internal state of homogenization scheme + crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity + forall ( i = 1:homogenization_maxNgrains, & + j = 1:mesh_maxNips, & + k = 1:mesh_NcpElems ) & + constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites + write(6,'(a10,/,4(3(f20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p + do k = 1,mesh_NcpElems + do j = 1,mesh_maxNips + if (homogenization_sizeState(j,k) > 0_pInt) & + homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme enddo enddo endif @@ -269,7 +263,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- case (3) - if (IP==1.AND.cp_en==1) write(6,*) 'Temp from CPFEM', Temperature materialpoint_Temperature(IP,cp_en) = Temperature materialpoint_F0(:,:,IP,cp_en) = ffn materialpoint_F(:,:,IP,cp_en) = ffn1 @@ -294,14 +287,13 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt ! return the local stress and the jacobian from storage cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en) jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en) + if (IP == 1 .and. cp_en == 1) write(6,'(a,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip 1 el 1',jacobian/1e9 ! return temperature if (theInc > 0_pInt) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. - return end subroutine - END MODULE CPFEM \ No newline at end of file diff --git a/code/IO.f90 b/code/IO.f90 index e5962c989..7e21af076 100644 --- a/code/IO.f90 +++ b/code/IO.f90 @@ -58,8 +58,9 @@ use prec, only: pReal, pInt implicit none + integer(pInt), intent(in) :: unit + integer(pInt) extPos character(256) outName - integer(pInt) unit, extPos character(3) ext inquire(6, name=outName) ! determine outputfileName @@ -78,6 +79,29 @@ endfunction +!******************************************************************** +! open (write) file related to current job +! but with different extension to given unit +!******************************************************************** + logical function IO_open_jobFile(unit,newExt) + + use prec, only: pReal, pInt + implicit none + + integer(pInt), intent(in) :: unit + character(*), intent(in) :: newExt + character(256) outName + + inquire(6, name=outName) ! determine outputfileName + open(unit,status='replace',err=100,file=outName(1:len_trim(outName)-3)//newExt) + IO_open_jobFile = .true. + return +100 IO_open_jobFile = .false. + return + + endfunction + + !******************************************************************** ! hybrid IA repetition counter !******************************************************************** @@ -771,6 +795,8 @@ endfunction select case (ID) case (0) msg = 'Unable to open input file' + case (50) + msg = 'Error writing constitutive output description' case (100) msg = 'Error reading from configuration file' case (105) @@ -793,22 +819,16 @@ endfunction msg = 'Unknown constitution specified' case (201) msg = 'Unknown homogenization specified' - case (202) - msg = 'Number of slip systems too small' - case (203) - msg = 'Negative initial slip resistance' - case (204) - msg = 'Non-positive reference shear rate' case (205) + msg = 'Unknown lattice structure encountered' + case (210) + msg = 'Negative initial resistance' + case (211) + msg = 'Non-positive reference shear rate' + case (212) msg = 'Non-positive stress exponent' - case (206) - msg = 'Non-positive initial hardening slope' - case (207) + case (213) msg = 'Non-positive saturation stress' - case (208) - msg = 'Non-positive w0' - case (209) - msg = 'Negative latent hardening ratio' case (220) msg = 'Negative initial dislocation density' case (221) @@ -818,9 +838,13 @@ endfunction case (223) msg = 'Negative self diffusion energy' case (224) - msg = 'Negative diffusion constant' + msg = 'Non-positive diffusion prefactor' + case (225) + msg = 'No slip systems specified' case (240) msg = 'Non-positive Taylor factor' + case (241) + msg = 'Non-positive hardening exponent' case (260) msg = 'Non-positive relevant strain' case (261) @@ -852,9 +876,7 @@ endfunction case (274) msg = 'Non-positive relative maximum value (upper bound) for GIA residual' case (275) - msg = 'Limit for GIA iteration too small' - case (276) - msg = 'Non-positive relative tolerance for temperature' + msg = 'Limit for GIA iteration too small' case (300) msg = 'This material can only be used with elements with three direct stress components' case (500) @@ -866,7 +888,7 @@ endfunction case (700) msg = 'Singular matrix in stress iteration' case (800) - msg = 'GIA requires 8 grains per IP (bonehead, you!)' + msg = 'RGC requires 8 grains per IP (bonehead, you!) -- but now outdated' case default msg = 'Unknown error number...' end select diff --git a/code/constitutive.f90 b/code/constitutive.f90 index b95708ca9..66b033a88 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -39,24 +39,64 @@ subroutine constitutive_init() !* Module initialization * !************************************** use prec, only: pReal,pInt - use IO, only: IO_error, IO_open_file + use debug, only: debugger + use IO, only: IO_error, IO_open_file, IO_open_jobFile use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use material - use constitutive_phenomenological use constitutive_j2 + use constitutive_phenopowerlaw use constitutive_dislobased integer(pInt), parameter :: fileunit = 200 - integer(pInt) e,i,g,myInstance,myNgrains - + integer(pInt) e,i,g,p,myInstance,myNgrains + integer(pInt), dimension(:,:), pointer :: thisSize + character(64), dimension(:,:), pointer :: thisOutput + logical knownConstitution + if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file - call constitutive_phenomenological_init(fileunit) ! parse all phases of this constitution - call constitutive_j2_init(fileunit) + call constitutive_j2_init(fileunit) ! parse all phases of this constitution + call constitutive_phenopowerlaw_init(fileunit) call constitutive_dislobased_init(fileunit) close(fileunit) +! write description file for constitutive phase output + + if(.not. IO_open_jobFile(fileunit,'outputConstitutive')) call IO_error (50) ! problems in writing file + + do p = 1,material_Nphase + i = phase_constitutionInstance(p) ! which instance of a constitution is present phase + knownConstitution = .true. ! assume valid + select case(phase_constitution(p)) ! split per constitiution + case (constitutive_j2_label) + thisOutput => constitutive_j2_output + thisSize => constitutive_j2_sizePostResult + case (constitutive_phenopowerlaw_label) + thisOutput => constitutive_phenopowerlaw_output + thisSize => constitutive_phenopowerlaw_sizePostResult + case (constitutive_dislobased_label) + thisOutput => constitutive_dislobased_output + thisSize => constitutive_dislobased_sizePostResult + case default + knownConstitution = .false. + end select + + write(fileunit,*) + write(fileunit,'(a)') '['//trim(phase_name(p))//']' + write(fileunit,*) + if (knownConstitution) then + write(fileunit,'(a)') '#'//char(9)//'constitution'//char(9)//trim(phase_constitution(p)) + do e = 1,phase_Noutput(p) + write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i) + enddo + endif + enddo + + close(fileunit) + +! allocate memory for state management + allocate(constitutive_state0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) @@ -69,17 +109,9 @@ subroutine constitutive_init() myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs do g = 1,myNgrains ! loop over grains + debugger = (e == 1 .and. i == 1 .and. g == 1) myInstance = phase_constitutionInstance(material_phase(g,i,e)) select case(phase_constitution(material_phase(g,i,e))) - case (constitutive_phenomenological_label) - allocate(constitutive_state0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) - allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) - allocate(constitutive_subState0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) - allocate(constitutive_state(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance))) - constitutive_state0(g,i,e)%p = constitutive_phenomenological_stateInit(myInstance) - constitutive_sizeDotState(g,i,e) = constitutive_phenomenological_sizeDotState(myInstance) - constitutive_sizeState(g,i,e) = constitutive_phenomenological_sizeState(myInstance) - constitutive_sizePostResults(g,i,e) = constitutive_phenomenological_sizePostResults(myInstance) case (constitutive_j2_label) allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) @@ -89,6 +121,15 @@ subroutine constitutive_init() constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) + case (constitutive_phenopowerlaw_label) + allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) + constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) + constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance) + constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) + constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance) case (constitutive_dislobased_label) allocate(constitutive_state0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) @@ -139,8 +180,8 @@ function constitutive_homogenizedC(ipc,ip,el) !********************************************************************* use prec, only: pReal,pInt use material, only: phase_constitution,material_phase - use constitutive_phenomenological use constitutive_j2 + use constitutive_phenopowerlaw use constitutive_dislobased implicit none @@ -149,13 +190,12 @@ function constitutive_homogenizedC(ipc,ip,el) real(pReal), dimension(6,6) :: constitutive_homogenizedC select case (phase_constitution(material_phase(ipc,ip,el))) - case (constitutive_phenomenological_label) - constitutive_homogenizedC = constitutive_phenomenological_homogenizedC(constitutive_state,ipc,ip,el) case (constitutive_j2_label) constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) + constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(constitutive_state,ipc,ip,el) case (constitutive_dislobased_label) constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state,ipc,ip,el) - end select return @@ -174,8 +214,8 @@ subroutine constitutive_microstructure(Temperature,ipc,ip,el) !********************************************************************* use prec, only: pReal,pInt use material, only: phase_constitution,material_phase - use constitutive_phenomenological use constitutive_j2 + use constitutive_phenopowerlaw use constitutive_dislobased implicit none @@ -184,13 +224,12 @@ integer(pInt) ipc,ip,el real(pReal) Temperature select case (phase_constitution(material_phase(ipc,ip,el))) - case (constitutive_phenomenological_label) - call constitutive_phenomenological_microstructure(Temperature,constitutive_state,ipc,ip,el) case (constitutive_j2_label) call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) + call constitutive_phenopowerlaw_microstructure(Temperature,constitutive_state,ipc,ip,el) case (constitutive_dislobased_label) call constitutive_dislobased_microstructure(Temperature,constitutive_state,ipc,ip,el) - end select end subroutine @@ -211,8 +250,8 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i !********************************************************************* use prec, only: pReal,pInt use material, only: phase_constitution,material_phase - use constitutive_phenomenological use constitutive_j2 + use constitutive_phenopowerlaw use constitutive_dislobased implicit none @@ -224,13 +263,12 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i real(pReal), dimension(9,9) :: dLp_dTstar select case (phase_constitution(material_phase(ipc,ip,el))) - case (constitutive_phenomenological_label) - call constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_j2_label) call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) + call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_dislobased_label) call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) - end select return @@ -252,8 +290,8 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el) !********************************************************************* use prec, only: pReal,pInt use material, only: phase_constitution,material_phase - use constitutive_phenomenological use constitutive_j2 + use constitutive_phenopowerlaw use constitutive_dislobased implicit none @@ -264,13 +302,12 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el) real(pReal), dimension(constitutive_sizeDotState(ipc,ip,el)) :: constitutive_dotState select case (phase_constitution(material_phase(ipc,ip,el))) - case (constitutive_phenomenological_label) - constitutive_dotState = constitutive_phenomenological_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_j2_label) constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) + constitutive_dotState = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_dislobased_label) constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) - end select return end function @@ -291,8 +328,8 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) !********************************************************************* use prec, only: pReal,pInt use material, only: phase_constitution,material_phase - use constitutive_phenomenological use constitutive_j2 + use constitutive_phenopowerlaw use constitutive_dislobased implicit none @@ -303,10 +340,10 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el) real(pReal) constitutive_dotTemperature select case (phase_constitution(material_phase(ipc,ip,el))) - case (constitutive_phenomenological_label) - constitutive_dotTemperature = constitutive_phenomenological_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_j2_label) constitutive_dotTemperature = constitutive_j2_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) + constitutive_dotTemperature = constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) case (constitutive_dislobased_label) constitutive_dotTemperature = constitutive_dislobased_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) @@ -327,8 +364,8 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el) !********************************************************************* use prec, only: pReal,pInt use material, only: phase_constitution,material_phase - use constitutive_phenomenological use constitutive_j2 + use constitutive_phenopowerlaw use constitutive_dislobased implicit none @@ -340,10 +377,10 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el) constitutive_postResults = 0.0_pReal select case (phase_constitution(material_phase(ipc,ip,el))) - case (constitutive_phenomenological_label) - constitutive_postResults = constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) case (constitutive_j2_label) constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) + case (constitutive_phenopowerlaw_label) + constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) case (constitutive_dislobased_label) constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) diff --git a/code/constitutive_dislobased.f90 b/code/constitutive_dislobased.f90 index 998e18453..c1f184578 100644 --- a/code/constitutive_dislobased.f90 +++ b/code/constitutive_dislobased.f90 @@ -8,48 +8,72 @@ !* - orientations * !************************************ +! [Alu] +! constitution dislobased +! (output) dislodensity +! (output) rateofshear +! lattice_structure 1 +! Nslip 12 +! +! c11 106.75e9 +! c12 60.41e9 +! c44 28.34e9 +! +! burgers 2.86e-10 # Burgers vector [m] +! Qedge 3e-19 # Activation energy for dislocation glide [J/K] (0.5*G*b^3) +! Qsd 2.4e-19 # Activation energy for self diffusion [J/K] (gamma-iron) +! diff0 1e-3 # prefactor vacancy diffusion coeffficent (gamma-iron) +! interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5 # Dislocation interaction coefficients +! +! rho0 6.0e12 # Initial dislocation density [m/m^3] +! +! c1 0.1 # Passing stress adjustment +! c2 2.0 # Jump width adjustment +! c3 1.0 # Activation volume adjustment +! c4 50.0 # Average slip distance adjustment for lock formation +! c7 8.0 # Athermal recovery adjustment +! c8 1.0e10 # Thermal recovery adjustment (plays no role for me) + MODULE constitutive_dislobased !*** Include other modules *** use prec, only: pReal,pInt implicit none character (len=*), parameter :: constitutive_dislobased_label = 'dislobased' - integer(pInt), dimension(:), allocatable :: constitutive_dislobased_sizeDotState, & - constitutive_dislobased_sizeState, & - constitutive_dislobased_sizePostResults - 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_Nslip - integer(pInt), dimension(:), allocatable :: constitutive_dislobased_Ntwin - real(pReal), dimension(:), allocatable :: constitutive_dislobased_C11 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_C12 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_C13 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_C33 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_C44 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Gmod - real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Cslip_66 - real(pReal), dimension(:,:,:,:), allocatable :: constitutive_dislobased_Ctwin_66 - real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_dislobased_Cslip_3333 - real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_dislobased_Ctwin_3333 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_rho0 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_bg - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Qedge - real(pReal), dimension(:), allocatable :: constitutive_dislobased_grainsize - real(pReal), dimension(:), allocatable :: constitutive_dislobased_stacksize - real(pReal), dimension(:), allocatable :: constitutive_dislobased_fmax - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Ndot0 - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Cmfpslip - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Cmfptwin - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Cthresholdslip - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Cthresholdtwin - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Cactivolume - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Cstorage - real(pReal), dimension(:), allocatable :: constitutive_dislobased_Carecovery - real(pReal), dimension(:), allocatable :: constitutive_dislobased_CoverA - real(pReal), dimension(:,:), allocatable :: constitutive_dislobased_SlipIntCoeff - real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Iparallel - real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Iforest + + integer(pInt), dimension(:), allocatable :: constitutive_dislobased_sizeDotState, & + constitutive_dislobased_sizeState, & + constitutive_dislobased_sizePostResults + integer(pInt), dimension(:,:), allocatable,target :: constitutive_dislobased_sizePostResult ! size of each post result output + character(len=64), dimension(:,:), allocatable,target :: constitutive_dislobased_output ! name of each post result output + character(len=32), dimension(:), allocatable :: constitutive_dislobased_structureName + integer(pInt), dimension(:), allocatable :: constitutive_dislobased_structure + integer(pInt), dimension(:), allocatable :: constitutive_dislobased_Nslip + real(pReal), dimension(:), allocatable :: constitutive_dislobased_C11 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_C12 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_C13 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_C33 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_C44 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_Gmod + real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Cslip_66 +!* Visco-plastic constitutive_phenomenological parameters + real(pReal), dimension(:), allocatable :: constitutive_dislobased_rho0 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_bg + real(pReal), dimension(:), allocatable :: constitutive_dislobased_Qedge + real(pReal), dimension(:), allocatable :: constitutive_dislobased_Qsd + real(pReal), dimension(:), allocatable :: constitutive_dislobased_D0 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_c1 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_c2 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_c3 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_c4 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_c5 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_c6 + real(pReal), dimension(:), allocatable :: constitutive_dislobased_c7 + 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_Iparallel + real(pReal), dimension(:,:,:), allocatable :: constitutive_dislobased_Iforest !************************************* !* Definition of material properties * @@ -70,7 +94,6 @@ CONTAINS !* - constitutive_microstructure !* - constitutive_LpAndItsTangent !* - consistutive_dotState -!* - constitutive_dotTemperature !* - consistutive_postResults !**************************************** @@ -82,50 +105,54 @@ subroutine constitutive_dislobased_init(file) use math, only: math_Mandel3333to66, math_Voigt66to3333, math_mul3x3 use IO use material - use lattice, only: lattice_sn,lattice_st,lattice_interactionSlipSlip,lattice_initializeStructure,lattice_Qtwin,lattice_tn + use lattice, only: lattice_sn, lattice_st, lattice_interactionSlipSlip, lattice_initializeStructure integer(pInt), intent(in) :: file integer(pInt), parameter :: maxNchunks = 7 integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) section,maxNinstance,i,j,k,l,m,n,o,p,q,r,s,output + integer(pInt) section, maxNinstance, i,j,k,l, output, mySize character(len=64) tag character(len=1024) line real(pReal) x,y + write(6,*) + write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_dislobased_label,' init -+>>>' + write(6,*) + maxNinstance = count(phase_constitution == constitutive_dislobased_label) if (maxNinstance == 0) return - allocate(constitutive_dislobased_sizeDotState(maxNinstance)) ; constitutive_dislobased_sizeDotState = 0_pInt - allocate(constitutive_dislobased_sizeState(maxNinstance)) ; constitutive_dislobased_sizeState = 0_pInt - allocate(constitutive_dislobased_sizePostResults(maxNinstance)) ; constitutive_dislobased_sizePostResults = 0_pInt - allocate(constitutive_dislobased_structureName(maxNinstance)) ; constitutive_dislobased_structureName = '' - allocate(constitutive_dislobased_structure(maxNinstance)) ; constitutive_dislobased_structure = 0_pInt - allocate(constitutive_dislobased_Nslip(maxNinstance)) ; constitutive_dislobased_Nslip = 0_pInt - allocate(constitutive_dislobased_Ntwin(maxNinstance)) ; constitutive_dislobased_Ntwin = 0_pInt - allocate(constitutive_dislobased_C11(maxNinstance)) ; constitutive_dislobased_C11 = 0.0_pReal - allocate(constitutive_dislobased_C12(maxNinstance)) ; constitutive_dislobased_C12 = 0.0_pReal - allocate(constitutive_dislobased_C13(maxNinstance)) ; constitutive_dislobased_C13 = 0.0_pReal - allocate(constitutive_dislobased_C33(maxNinstance)) ; constitutive_dislobased_C33 = 0.0_pReal - allocate(constitutive_dislobased_C44(maxNinstance)) ; constitutive_dislobased_C44 = 0.0_pReal - allocate(constitutive_dislobased_Gmod(maxNinstance)) ; constitutive_dislobased_Gmod = 0.0_pReal - allocate(constitutive_dislobased_Cslip_66(6,6,maxNinstance)) ; constitutive_dislobased_Cslip_66 = 0.0_pReal - allocate(constitutive_dislobased_Cslip_3333(3,3,3,3,maxNinstance)) ; constitutive_dislobased_Ctwin_3333 = 0.0_pReal - allocate(constitutive_dislobased_rho0(maxNinstance)) ; constitutive_dislobased_rho0 = 0.0_pReal - allocate(constitutive_dislobased_bg(maxNinstance)) ; constitutive_dislobased_bg = 0.0_pReal - allocate(constitutive_dislobased_Qedge(maxNinstance)) ; constitutive_dislobased_Qedge = 0.0_pReal - allocate(constitutive_dislobased_grainsize(maxNinstance)) ; constitutive_dislobased_grainsize = 0.0_pReal - allocate(constitutive_dislobased_stacksize(maxNinstance)) ; constitutive_dislobased_stacksize = 0.0_pReal - allocate(constitutive_dislobased_fmax(maxNinstance)) ; constitutive_dislobased_fmax = 0.0_pReal - allocate(constitutive_dislobased_Ndot0(maxNinstance)) ; constitutive_dislobased_Ndot0 = 0.0_pReal - allocate(constitutive_dislobased_Cmfpslip(maxNinstance)) ; constitutive_dislobased_Cmfpslip = 0.0_pReal - allocate(constitutive_dislobased_Cmfptwin(maxNinstance)) ; constitutive_dislobased_Cmfptwin = 0.0_pReal - allocate(constitutive_dislobased_Cthresholdslip(maxNinstance)) ; constitutive_dislobased_Cthresholdslip = 0.0_pReal - allocate(constitutive_dislobased_Cthresholdtwin(maxNinstance)) ; constitutive_dislobased_Cthresholdtwin = 0.0_pReal - allocate(constitutive_dislobased_Cactivolume(maxNinstance)) ; constitutive_dislobased_Cactivolume = 0.0_pReal - allocate(constitutive_dislobased_Cstorage(maxNinstance)) ; constitutive_dislobased_Cstorage = 0.0_pReal - allocate(constitutive_dislobased_Carecovery(maxNinstance)) ; constitutive_dislobased_Carecovery = 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_output(maxval(phase_Noutput),maxNinstance)) ; constitutive_dislobased_output = '' + allocate(constitutive_dislobased_sizeDotState(maxNinstance)) ; constitutive_dislobased_sizeDotState = 0_pInt + allocate(constitutive_dislobased_sizeState(maxNinstance)) ; constitutive_dislobased_sizeState = 0_pInt + allocate(constitutive_dislobased_sizePostResults(maxNinstance)); constitutive_dislobased_sizePostResults = 0_pInt + allocate(constitutive_dislobased_sizePostResult(maxval(phase_Noutput), & + maxNinstance)) ; constitutive_dislobased_sizePostResult = 0_pInt + allocate(constitutive_dislobased_output(maxval(phase_Noutput), & + 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_Nslip(maxNinstance)) ; constitutive_dislobased_Nslip = 0_pInt + allocate(constitutive_dislobased_C11(maxNinstance)) ; constitutive_dislobased_C11 = 0.0_pReal + allocate(constitutive_dislobased_C12(maxNinstance)) ; constitutive_dislobased_C12 = 0.0_pReal + allocate(constitutive_dislobased_C13(maxNinstance)) ; constitutive_dislobased_C13 = 0.0_pReal + allocate(constitutive_dislobased_C33(maxNinstance)) ; constitutive_dislobased_C33 = 0.0_pReal + allocate(constitutive_dislobased_C44(maxNinstance)) ; constitutive_dislobased_C44 = 0.0_pReal + allocate(constitutive_dislobased_Gmod(maxNinstance)) ; constitutive_dislobased_Gmod = 0.0_pReal + allocate(constitutive_dislobased_Cslip_66(6,6,maxNinstance)) ; constitutive_dislobased_Cslip_66 = 0.0_pReal + allocate(constitutive_dislobased_rho0(maxNinstance)) ; constitutive_dislobased_rho0 = 0.0_pReal + allocate(constitutive_dislobased_bg(maxNinstance)) ; constitutive_dislobased_bg = 0.0_pReal + allocate(constitutive_dislobased_Qedge(maxNinstance)) ; constitutive_dislobased_Qedge = 0.0_pReal + allocate(constitutive_dislobased_Qsd(maxNinstance)) ; constitutive_dislobased_Qsd = 0.0_pReal + allocate(constitutive_dislobased_D0(maxNinstance)) ; constitutive_dislobased_D0 = 0.0_pReal + allocate(constitutive_dislobased_c1(maxNinstance)) ; constitutive_dislobased_c1 = 0.0_pReal + allocate(constitutive_dislobased_c2(maxNinstance)) ; constitutive_dislobased_c2 = 0.0_pReal + allocate(constitutive_dislobased_c3(maxNinstance)) ; constitutive_dislobased_c3 = 0.0_pReal + allocate(constitutive_dislobased_c4(maxNinstance)) ; constitutive_dislobased_c4 = 0.0_pReal + allocate(constitutive_dislobased_c5(maxNinstance)) ; constitutive_dislobased_c5 = 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_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 rewind(file) line = '' @@ -157,8 +184,6 @@ subroutine constitutive_dislobased_init(file) constitutive_dislobased_CoverA(i) = IO_floatValue(line,positions,2) case ('nslip') constitutive_dislobased_Nslip(i) = IO_intValue(line,positions,2) - case ('ntwin') - constitutive_dislobased_Ntwin(i) = IO_intValue(line,positions,2) case ('c11') constitutive_dislobased_C11(i) = IO_floatValue(line,positions,2) case ('c12') @@ -175,31 +200,29 @@ subroutine constitutive_dislobased_init(file) constitutive_dislobased_bg(i) = IO_floatValue(line,positions,2) case ('qedge') constitutive_dislobased_Qedge(i) = IO_floatValue(line,positions,2) - case ('grainsize') - constitutive_dislobased_grainsize(i) = IO_floatValue(line,positions,2) - case ('stacksize') - constitutive_dislobased_stacksize(i) = IO_floatValue(line,positions,2) - case ('fmax') - constitutive_dislobased_fmax(i) = IO_floatValue(line,positions,2) - case ('ndot0') - constitutive_dislobased_Ndot0(i) = IO_floatValue(line,positions,2) - case ('cmfpslip') - constitutive_dislobased_Cmfpslip(i) = IO_floatValue(line,positions,2) - case ('cmfptwin') - constitutive_dislobased_Cmfptwin(i) = IO_floatValue(line,positions,2) - case ('cthresholdslip') - constitutive_dislobased_Cthresholdslip(i) = IO_floatValue(line,positions,2) - case ('cthresholdtwin') - constitutive_dislobased_Cthresholdtwin(i) = IO_floatValue(line,positions,2) - case ('cactivolume') - constitutive_dislobased_Cactivolume(i) = IO_floatValue(line,positions,2) - case ('cstorage') - constitutive_dislobased_Cstorage(i) = IO_floatValue(line,positions,2) - case ('carecovery') - constitutive_dislobased_Carecovery(i) = IO_floatValue(line,positions,2) + case ('qsd') + constitutive_dislobased_Qsd(i) = IO_floatValue(line,positions,2) + case ('diff0') + constitutive_dislobased_D0(i) = IO_floatValue(line,positions,2) + case ('c1') + constitutive_dislobased_c1(i) = IO_floatValue(line,positions,2) + case ('c2') + constitutive_dislobased_c2(i) = IO_floatValue(line,positions,2) + case ('c3') + constitutive_dislobased_c3(i) = IO_floatValue(line,positions,2) + case ('c4') + constitutive_dislobased_c4(i) = IO_floatValue(line,positions,2) + case ('c5') + constitutive_dislobased_c5(i) = IO_floatValue(line,positions,2) + case ('c6') + constitutive_dislobased_c6(i) = IO_floatValue(line,positions,2) + case ('c7') + constitutive_dislobased_c7(i) = IO_floatValue(line,positions,2) + case ('c8') + constitutive_dislobased_c8(i) = IO_floatValue(line,positions,2) case ('interaction_coefficients') - forall (j=2:min(7,positions(1))) & - constitutive_dislobased_SlipIntCoeff(j-1,i) = IO_floatValue(line,positions,j) + forall (j=1:6) & + constitutive_dislobased_SlipIntCoeff(j,i) = IO_floatValue(line,positions,1+j) end select endif enddo @@ -209,46 +232,43 @@ subroutine constitutive_dislobased_init(file) constitutive_dislobased_structure(i) = lattice_initializeStructure(constitutive_dislobased_structureName(i), & 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_structure(i) < 1) call IO_error(205) 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_Qedge(i) <= 0.0_pReal) call IO_error(222) + if (constitutive_dislobased_Qsd(i) <= 0.0_pReal) call IO_error(223) + if (constitutive_dislobased_D0(i) <= 0.0_pReal) call IO_error(224) + if (constitutive_dislobased_Nslip(i) < 1) call IO_error(225) enddo - allocate(constitutive_dislobased_Iparallel(maxval(constitutive_dislobased_Nslip),maxval(constitutive_dislobased_Nslip),maxNinstance)) - constitutive_dislobased_Iparallel = 0.0_pReal - allocate(constitutive_dislobased_Iforest(maxval(constitutive_dislobased_Nslip),maxval(constitutive_dislobased_Nslip),maxNinstance)) - constitutive_dislobased_Iforest = 0.0_pReal - allocate(constitutive_dislobased_Ctwin_66(6,6,maxval(constitutive_dislobased_Ntwin),maxNinstance)) - constitutive_dislobased_Ctwin_66 = 0.0_pReal - allocate(constitutive_dislobased_Ctwin_3333(3,3,3,3,maxval(constitutive_dislobased_Ntwin),maxNinstance)) - constitutive_dislobased_Ctwin_3333 = 0.0_pReal + allocate(constitutive_dislobased_Iparallel(maxval(constitutive_dislobased_Nslip),& + maxval(constitutive_dislobased_Nslip),& + maxNinstance)) - do i = 1,maxNinstance - constitutive_dislobased_sizeDotState(i) = constitutive_dislobased_Nslip(i)+constitutive_dislobased_Ntwin(i) - constitutive_dislobased_sizeState(i) = 10*constitutive_dislobased_Nslip(i)+5*constitutive_dislobased_Ntwin(i) + allocate(constitutive_dislobased_Iforest(maxval(constitutive_dislobased_Nslip),& + maxval(constitutive_dislobased_Nslip),& + maxNinstance)) + do i = 1,maxNinstance do j = 1,maxval(phase_Noutput) - select case(constitutive_dislobased_output(j,i)) - case('state_slip') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Nslip(i) - case('rateofshear_slip') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Nslip(i) - case('mfp_slip') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Nslip(i) - case('thresholdstress_slip') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Nslip(i) - case('state_twin') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Ntwin(i) - case('rateofshear_twin') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Ntwin(i) - case('mfp_twin') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Ntwin(i) - case('thresholdstress_twin') - constitutive_dislobased_sizePostResults(i) = constitutive_dislobased_sizePostResults(i) + constitutive_dislobased_Ntwin(i) - end select + select case(constitutive_dislobased_output(j,i)) + case('dislodensity') + mySize = constitutive_dislobased_Nslip(i) + case('rateofshear') + mySize = constitutive_dislobased_Nslip(i) + case default + mySize = 0_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + constitutive_dislobased_sizePostResult(j,i) = mySize + constitutive_dislobased_sizePostResults(i) = & + constitutive_dislobased_sizePostResults(i) + mySize + endif enddo + + constitutive_dislobased_sizeDotState(i) = constitutive_dislobased_Nslip(i) + constitutive_dislobased_sizeState(i) = 8*constitutive_dislobased_Nslip(i) constitutive_dislobased_Gmod(i) = constitutive_dislobased_C44(i) select case (constitutive_dislobased_structure(i)) @@ -274,33 +294,9 @@ subroutine constitutive_dislobased_init(file) constitutive_dislobased_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_dislobased_C11(i)- & constitutive_dislobased_C12(i)) end select - constitutive_dislobased_Cslip_66(:,:,i) = math_Mandel3333to66(math_Voigt66to3333(constitutive_dislobased_Cslip_66(:,:,i))) - - !* Construction of the twin elasticity matrices - !* Iteration over the systems - constitutive_dislobased_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_dislobased_Cslip_66(:,:,i)) - do j=1,constitutive_dislobased_Ntwin(i) - do k=1,3 - do l=1,3 - do m=1,3 - do n=1,3 - do p=1,3 - do q=1,3 - do r=1,3 - do s=1,3 - constitutive_dislobased_Ctwin_3333(k,l,m,n,j,i) = constitutive_dislobased_Ctwin_3333(k,l,m,n,j,i) + & - constitutive_dislobased_Cslip_3333(p,q,r,s,i)*& - lattice_Qtwin(k,p,j,i)*lattice_Qtwin(l,q,j,i)*lattice_Qtwin(m,r,j,i)*lattice_Qtwin(n,s,j,i) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - enddo - constitutive_dislobased_Ctwin_66(:,:,j,i) = math_Mandel3333to66(constitutive_dislobased_Ctwin_3333(:,:,:,:,j,i)) - enddo + constitutive_dislobased_Cslip_66(:,:,i) = & + math_Mandel3333to66(math_Voigt66to3333(constitutive_dislobased_Cslip_66(:,:,i))) + !* Construction of the hardening matrices !* Iteration over the systems @@ -321,6 +317,7 @@ subroutine constitutive_dislobased_init(file) enddo return + end subroutine @@ -331,193 +328,108 @@ function constitutive_dislobased_stateInit(myInstance) use prec, only: pReal,pInt implicit none - !* Definition of variables +!* Definition of variables integer(pInt), intent(in) :: myInstance - real(pReal), dimension(constitutive_dislobased_sizeState(myInstance)) :: constitutive_dislobased_stateInit - - constitutive_dislobased_stateInit = 0.0_pReal - constitutive_dislobased_stateInit(1:constitutive_dislobased_Nslip(myInstance)) = constitutive_dislobased_rho0(myInstance) + real(pReal), dimension(constitutive_dislobased_Nslip(myInstance)) :: constitutive_dislobased_stateInit + + constitutive_dislobased_stateInit = constitutive_dislobased_rho0(myInstance) return end function - function constitutive_dislobased_homogenizedC(state,ipc,ip,el) !********************************************************************* -!* calculates homogenized elacticity matrix * -!* - state : microstructure quantities * +!* homogenized elacticity matrix * +!* INPUT: * +!* - state : state variables * !* - ipc : component-ID of current integration point * !* - ip : current integration point * !* - el : current element * !********************************************************************* - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_NcpElems,mesh_maxNips + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance implicit none - !* Definition of variables +!* Definition of variables integer(pInt), intent(in) :: ipc,ip,el - integer(pInt) matID,ns,nt,i - real(pReal) sumf + integer(pInt) matID real(pReal), dimension(6,6) :: constitutive_dislobased_homogenizedC type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - !* Shortened notation matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - ns = constitutive_dislobased_Nslip(matID) - nt = constitutive_dislobased_Ntwin(matID) - - !* Total twin volume fraction - sumf = 0.0_pReal - if (nt > 0_pInt) sumf = sum(state(ipc,ip,el)%p((ns+1):(ns+nt))) - - !* Homogenized elasticity matrix - constitutive_dislobased_homogenizedC = (1.0_pReal-sumf)*constitutive_dislobased_Cslip_66(:,:,matID) - do i=1,nt - constitutive_dislobased_homogenizedC = constitutive_dislobased_homogenizedC + & - state(ipc,ip,el)%p(ns+i)*constitutive_dislobased_Ctwin_66(:,:,i,matID) - enddo + constitutive_dislobased_homogenizedC = constitutive_dislobased_Cslip_66(:,:,matID) return + end function subroutine constitutive_dislobased_microstructure(Temperature,state,ipc,ip,el) !********************************************************************* -!* calculates quantities characterizing the microstructure * -!* - Temperature : temperature * -!* - state : microstructure quantities * +!* calculate derived quantities from state (not used here) * +!* INPUT: * +!* - Tp : temperature * !* - ipc : component-ID of current integration point * !* - ip : current integration point * !* - el : current element * !********************************************************************* - use prec, only: pReal,pInt,p_vec - use math, only: pi - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance - use lattice, only: lattice_interactionSlipTwin,lattice_interactionTwinTwin + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance implicit none - !* Definition of variables - integer(pInt), intent(in) :: ipc,ip,el - integer(pInt) matID,ns,nt,i - real(pReal) Temperature,sumf +!* Definition of variables + integer(pInt) ipc,ip,el,matID,n,i + real(pReal) Temperature type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - !* Shortened notation matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - ns = constitutive_dislobased_Nslip(matID) - nt = constitutive_dislobased_Ntwin(matID) - !* State: 1 : ns rho_ssd - !* State: ns+1 : ns+nt f - !* State: ns+nt+1 : 2*ns+nt rho_forest - !* State: 2*ns+nt+1 : 3*ns+nt rho_parallel - !* State: 3*ns+nt+1 : 4*ns+nt 1/lambda_slip - !* State: 4*ns+nt+1 : 5*ns+nt 1/lambda_sliptwin - !* State: 5*ns+nt+1 : 5*ns+2*nt 1/lambda_twin - !* State: 5*ns+2*nt+1 : 6*ns+2*nt mfp_slip - !* State: 6*ns+2*nt+1 : 6*ns+3*nt mfp_twin - !* State: 6*ns+3*nt+1 : 7*ns+3*nt threshold_stress_slip - !* State: 7*ns+3*nt+1 : 7*ns+4*nt threshold_stress_twin - !* State: 7*ns+4*nt+1 : 8*ns+4*nt activation volume - !* State: 8*ns+4*nt+1 : 8*ns+5*nt twin volume - !* State: 8*ns+5*nt+1 : 9*ns+5*nt rho_mobile - !* State: 9*ns+5*nt+1 : 10*ns+5*nt initial shear rate - - !* Total twin volume fraction - sumf = 0.0_pReal - if (nt > 0_pInt) sumf = sum(state(ipc,ip,el)%p((ns+1):(ns+nt))) - - !* Forest and parallel dislocation densities + n = constitutive_dislobased_Nslip(matID) + !* Quantities derived from state - slip + !* State: 1 : n rho + !* n+1 : 2n rho_f + !* 2n+1 : 3n rho_p + !* 3n+1 : 4n passing stress + !* 4n+1 : 5n jump width + !* 5n+1 : 6n activation volume + !* 6n+1 : 7n rho_m + !* 7n+1 : 8n g0_slip !$OMP CRITICAL (evilmatmul) - state(ipc,ip,el)%p((ns+nt+1):(2*ns+nt)) = matmul(constitutive_dislobased_Iforest (1:ns,1:ns,matID),state(ipc,ip,el)%p(1:ns)) - state(ipc,ip,el)%p((2*ns+nt+1):(3*ns+nt)) = matmul(constitutive_dislobased_Iparallel(1:ns,1:ns,matID),state(ipc,ip,el)%p(1:ns)) - !$OMP END CRITICAL (evilmatmul) - - !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation - do i=1,ns - state(ipc,ip,el)%p(3*ns+nt+i) = sqrt(state(ipc,ip,el)%p(ns+nt+i)) - enddo - - !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation - !$OMP CRITICAL (evilmatmul) - state(ipc,ip,el)%p((4*ns+nt+1):(5*ns+nt)) = 0.0_pReal - if (nt > 0_pInt) state(ipc,ip,el)%p((4*ns+nt+1):(5*ns+nt)) = & - matmul(lattice_interactionSlipTwin(1:ns,1:nt,constitutive_dislobased_structure(matID)),state(ipc,ip,el)%p((ns+1):(ns+nt)))/& - (2.0_pReal*constitutive_dislobased_stacksize(matID)*(1.0_pReal-sumf)) + state(ipc,ip,el)%p((n+1):(2*n)) = matmul(constitutive_dislobased_Iforest (1:n,1:n,matID),state(ipc,ip,el)%p(1:n)) + state(ipc,ip,el)%p((2*n+1):(3*n)) = matmul(constitutive_dislobased_Iparallel(1:n,1:n,matID),state(ipc,ip,el)%p(1:n)) !$OMP END CRITICAL (evilmatmul) - !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin - !$OMP CRITICAL (evilmatmul) - if (nt > 0_pInt) state(ipc,ip,el)%p((5*ns+nt+1):(5*ns+2*nt)) = & - matmul(lattice_interactionTwinTwin(1:nt,1:nt,constitutive_dislobased_structure(matID)),state(ipc,ip,el)%p((ns+1):(ns+nt)))/& - (2.0_pReal*constitutive_dislobased_stacksize(matID)*(1.0_pReal-sumf)) - !$OMP END CRITICAL (evilmatmul) + do i=1,n + + state(ipc,ip,el)%p(3*n+i) = & + constitutive_dislobased_c1(matID)*constitutive_dislobased_Gmod(matID)*& + constitutive_dislobased_bg(matID)*sqrt(state(ipc,ip,el)%p(2*n+i)) + + state(ipc,ip,el)%p(4*n+i) = & + constitutive_dislobased_c2(matID)/sqrt(state(ipc,ip,el)%p(n+i)) + + state(ipc,ip,el)%p(5*n+i) = & + constitutive_dislobased_c3(matID)*state(ipc,ip,el)%p(4*n+i)*constitutive_dislobased_bg(matID)**2.0_pReal + + state(ipc,ip,el)%p(6*n+i) = & + (2.0_pReal*kB*Temperature*sqrt(state(ipc,ip,el)%p(2*n+i)))/& + (constitutive_dislobased_c1(matID)*constitutive_dislobased_c3(matID)*constitutive_dislobased_Gmod(matID)*& + state(ipc,ip,el)%p(4*n+i)*constitutive_dislobased_bg(matID)**3.0_pReal) + + 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)*& + exp(-constitutive_dislobased_Qedge(matID)/(kB*Temperature)) - !* mean free path between 2 obstacles seen by a moving dislocation - do i=1,ns - if (nt > 0_pInt) then - state(ipc,ip,el)%p(5*ns+2*nt+i) = (constitutive_dislobased_Cmfpslip(matID)*constitutive_dislobased_grainsize(matID))/& - (1.0_pReal+constitutive_dislobased_grainsize(matID)*& - (state(ipc,ip,el)%p(3*ns+nt+i)+state(ipc,ip,el)%p(4*ns+nt+i))) - else - state(ipc,ip,el)%p(5*ns+i) = (constitutive_dislobased_Cmfpslip(matID)*constitutive_dislobased_grainsize(matID))/& - (1.0_pReal+constitutive_dislobased_grainsize(matID)*(state(ipc,ip,el)%p(3*ns+i))) - endif enddo - !* mean free path between 2 obstacles seen by a growing twin - do i=1,nt - state(ipc,ip,el)%p(6*ns+2*nt+i) = (constitutive_dislobased_Cmfptwin(matID)*constitutive_dislobased_grainsize(matID))/& - (1.0_pReal+constitutive_dislobased_grainsize(matID)*state(ipc,ip,el)%p(5*ns+nt+i)) - enddo - - !* threshold stress for dislocation motion - do i=1,ns - state(ipc,ip,el)%p(6*ns+3*nt+i) = constitutive_dislobased_Cthresholdslip(matID)*& - constitutive_dislobased_bg(matID)*constitutive_dislobased_Gmod(matID)*sqrt(state(ipc,ip,el)%p(2*ns+nt+i)) - enddo - - !* threshold stress for growing twin - do i=1,nt - state(ipc,ip,el)%p(7*ns+3*nt+i) = constitutive_dislobased_Cthresholdtwin(matID)*(sqrt(3.0_pReal)/3.0_pReal)*(& - (0.0002_pReal*Temperature-0.0396_pReal)/constitutive_dislobased_bg(matID)+& - (constitutive_dislobased_bg(matID)*constitutive_dislobased_Gmod(matID))/state(ipc,ip,el)%p(5*ns+2*nt+i)) - enddo - - !* activation volume for dislocation glide - do i=1,ns - state(ipc,ip,el)%p(7*ns+4*nt+i) = constitutive_dislobased_Cactivolume(matID)*& - constitutive_dislobased_bg(matID)*constitutive_dislobased_bg(matID)*state(ipc,ip,el)%p(5*ns+2*nt+i) - enddo - - !* final twin volume after growth - do i=1,nt - state(ipc,ip,el)%p(8*ns+4*nt+i) = (pi/6.0_pReal)*constitutive_dislobased_stacksize(matID)*& - state(ipc,ip,el)%p(6*ns+2*nt+i)*state(ipc,ip,el)%p(6*ns+2*nt+i) - enddo - - !* mobile dislocation densities - do i=1,ns - state(ipc,ip,el)%p(8*ns+5*nt+i) = (2.0_pReal*kB*Temperature*state(ipc,ip,el)%p(2*ns+nt+i))/& - (state(ipc,ip,el)%p(6*ns+3*nt+i)*state(ipc,ip,el)%p(7*ns+4*nt+i)) - enddo - - !* initial shear rate for slip - do i=1,ns - state(ipc,ip,el)%p(9*ns+5*nt+i) = state(ipc,ip,el)%p(8*ns+5*nt+i)*constitutive_dislobased_bg(matID)*attack_frequency*& - state(ipc,ip,el)%p(5*ns+2*nt+i)*exp(-constitutive_dislobased_Qedge(matID)/(kB*Temperature)) - enddo - end subroutine subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el) !********************************************************************* -!* calculates plastic velocity gradient and its tangent * +!* plastic velocity gradient and its tangent * !* INPUT: * -!* - Temperature : temperature * -!* - state : microstructure quantities * !* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * !* - ipc : component-ID at current integration point * !* - ip : current integration point * @@ -526,17 +438,18 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera !* - Lp : plastic velocity gradient * !* - dLp_dTstar : derivative of Lp (4th-rank tensor) * !********************************************************************* - use prec, only: pReal,pInt,p_vec - use math, only: math_Plain3333to99 - use mesh, only: mesh_NcpElems,mesh_maxNips + use prec, only: pReal,pInt,p_vec + use math, only: math_Plain3333to99, math_mul6x6 + use lattice, only: lattice_Sslip,lattice_Sslip_v + use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance - use lattice, only: lattice_Sslip,lattice_Stwin,lattice_Sslip_v,lattice_Stwin_v,lattice_shearTwin + implicit none - !* Definition of variables +!* Definition of variables integer(pInt) ipc,ip,el - integer(pInt) matID,i,k,l,m,n,ns,nt - real(pReal) Temperature,sumf + integer(pInt) matID,i,k,l,m,n + real(pReal) Temperature type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(3,3) :: Lp @@ -544,71 +457,34 @@ subroutine constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Tempera real(pReal), dimension(9,9) :: dLp_dTstar real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & gdot_slip,dgdot_dtauslip,tau_slip - real(pReal), dimension(constitutive_dislobased_Ntwin(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & - gdot_twin,dgdot_dtautwin,tau_twin - !* Shortened notation matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - ns = constitutive_dislobased_Nslip(matID) - nt = constitutive_dislobased_Ntwin(matID) + n = constitutive_dislobased_Nslip(matID) - !* Total twin volume fraction - sumf = 0.0_pReal - if (nt > 0_pInt) sumf = sum(state(ipc,ip,el)%p((ns+1):(ns+nt))) - - !* Calculation of Lp from dislocation glide +!* Calculation of Lp Lp = 0.0_pReal gdot_slip = 0.0_pReal - do i = 1,ns - tau_slip(i) = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) - - if ( abs(tau_slip(i)) > state(ipc,ip,el)%p(6*ns+3*nt+i) ) & - gdot_slip(i) = state(ipc,ip,el)%p(9*ns+5*nt+i)*sign(1.0_pReal,tau_slip(i))*& - sinh(((abs(tau_slip(i))-state(ipc,ip,el)%p(6*ns+3*nt+i))*state(ipc,ip,el)%p(7*ns+4*nt+i))/(kB*Temperature)) + do i = 1,constitutive_dislobased_Nslip(matID) + tau_slip(i) = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) + if ((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))>0) & + 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 + (1.0_pReal - sumf)*gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_dislobased_structure(matID)) + Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_dislobased_structure(matID)) enddo - !* Calculation of Lp from deformation twinning - gdot_twin = 0.0_pReal - do i = 1,nt - tau_twin(i) = dot_product(Tstar_v,lattice_Stwin_v(:,i,constitutive_dislobased_structure(matID))) - - if ( tau_twin(i) > 0.0_pReal ) & - gdot_twin(i) = (constitutive_dislobased_fmax(matID) - sumf)*lattice_shearTwin(i,constitutive_dislobased_structure(matID))*& - state(ipc,ip,el)%p(8*ns+4*nt+i)*constitutive_dislobased_Ndot0(matID)*& - exp(-(state(ipc,ip,el)%p(7*ns+3*nt+i)/tau_twin(i))**10.0_pReal) - - Lp = Lp + gdot_twin(i)*lattice_Stwin(:,:,i,constitutive_dislobased_structure(matID)) - enddo - - !* Calculation of the tangent of Lp from dislocation glide +!* Calculation of the tangent of Lp dLp_dTstar3333 = 0.0_pReal dLp_dTstar = 0.0_pReal dgdot_dtauslip = 0.0_pReal - do i = 1,ns - - if ( abs(tau_slip(i)) > state(ipc,ip,el)%p(6*ns+3*nt+i) ) & - dgdot_dtauslip(i) = (state(ipc,ip,el)%p(9*ns+5*nt+i)*state(ipc,ip,el)%p(7*ns+4*nt+i))/(kB*Temperature)*& - cosh(((abs(tau_slip(i))-state(ipc,ip,el)%p(6*ns+3*nt+i))*state(ipc,ip,el)%p(7*ns+4*nt+i))/(kB*Temperature)) - - 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) + & - dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_dislobased_structure(matID))& - *lattice_Sslip(m,n,i,constitutive_dislobased_structure(matID)) - enddo - - !* Calculation of the tangent of Lp from deformation twinning - dgdot_dtautwin = 0.0_pReal - do i = 1,nt - - if ( tau_twin(i) > 0.0_pReal ) & - dgdot_dtautwin(i) = (gdot_twin(i)*10.0_pReal*state(ipc,ip,el)%p(7*ns+3*nt+i)**10.0_pReal)/(tau_twin(i)**11.0_pReal) - - 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) + & - dgdot_dtautwin(i)*lattice_Stwin(k,l,i,constitutive_dislobased_structure(matID)) & - *lattice_Stwin(m,n,i,constitutive_dislobased_structure(matID)) + do i = 1,constitutive_dislobased_Nslip(matID) + if ((abs(tau_slip(i))-state(ipc,ip,el)%p(3*n+i))>0) & + 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) & + dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & + dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_dislobased_structure(matID))* & + lattice_Sslip(m,n,i,constitutive_dislobased_structure(matID)) enddo dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) @@ -620,8 +496,6 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el) !********************************************************************* !* rate of change of microstructure * !* INPUT: * -!* - Temperature : temperature * -!* - state : microstructure quantities * !* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * !* - ipc : component-ID at current integration point * !* - ip : current integration point * @@ -629,96 +503,80 @@ function constitutive_dislobased_dotState(Tstar_v,Temperature,state,ipc,ip,el) !* OUTPUT: * !* - constitutive_dotState : evolution of state variable * !********************************************************************* - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance - use lattice, only: lattice_Sslip_v,lattice_Stwin_v - implicit none - -!* Definition of variables - integer(pInt) ipc,ip,el - integer(pInt) matID,i,ns,nt - real(pReal) Temperature,sumf,tau_slip,tau_twin,gdot_slip,gdot_twin,storage,arecovery - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - real(pReal), dimension(6) :: Tstar_v - real(pReal), dimension(constitutive_dislobased_sizeDotState(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & - constitutive_dislobased_dotState - - !* Shortened notation - matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - ns = constitutive_dislobased_Nslip(matID) - nt = constitutive_dislobased_Ntwin(matID) - - !* Total twin volume fraction - sumf = 0.0_pReal - if (nt > 0_pInt) sumf = sum(state(ipc,ip,el)%p((ns+1):(ns+nt))) - - !* Dislocation density evolution - constitutive_dislobased_dotState = 0.0_pReal - do i = 1,ns - - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) - - if ( abs(tau_slip) > state(ipc,ip,el)%p(6*ns+3*nt+i) ) then - gdot_slip = state(ipc,ip,el)%p(9*ns+5*nt+i)*sign(1.0_pReal,tau_slip)*& - sinh(((abs(tau_slip)-state(ipc,ip,el)%p(6*ns+3*nt+i))*state(ipc,ip,el)%p(7*ns+4*nt+i))/(kB*Temperature) ) - - storage = (constitutive_dislobased_Cstorage(matID)*abs(gdot_slip))/& - (constitutive_dislobased_bg(matID)*state(ipc,ip,el)%p(5*ns+2*nt+i)) - - arecovery = constitutive_dislobased_Carecovery(matID)*state(ipc,ip,el)%p(i)*abs(gdot_slip) - - constitutive_dislobased_dotState(i) = storage - arecovery - endif - enddo - - !* Twin volume fraction evolution - do i = 1,nt - - tau_twin = dot_product(Tstar_v,lattice_Stwin_v(:,i,constitutive_dislobased_structure(matID))) - - if ( tau_twin > 0.0_pReal ) & - constitutive_dislobased_dotState(ns+i) = (constitutive_dislobased_fmax(matID) - sumf)*& - state(ipc,ip,el)%p(8*ns+4*nt+i)*constitutive_dislobased_Ndot0(matID)*& - exp(-(state(ipc,ip,el)%p(7*ns+3*nt+i)/tau_twin)**10.0_pReal) - - enddo - - return -end function - - -function constitutive_dislobased_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el) -!********************************************************************* -!* rate of change of microstructure * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - ipc : component-ID at current integration point * -!* - ip : current integration point * -!* - el : current element * -!* OUTPUT: * -!* - constitutive_dotTemperature : evolution of Temperature * -!********************************************************************* - use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v - use mesh, only: mesh_NcpElems,mesh_maxNips + use prec, only: pReal,pInt,p_vec + use lattice, only: lattice_Sslip_v + use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance implicit none !* Definition of variables integer(pInt) ipc,ip,el integer(pInt) matID,i,n - real(pReal) Temperature + real(pReal) Temperature,tau_slip,gdot_slip,locks,athermal_recovery,thermal_recovery type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state real(pReal), dimension(6) :: Tstar_v - real(pReal) constitutive_dislobased_dotTemperature - - constitutive_dislobased_dotTemperature = 0.0_pReal - + real(pReal), dimension(constitutive_dislobased_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + constitutive_dislobased_dotState + + matID = phase_constitutionInstance(material_phase(ipc,ip,el)) + n = constitutive_dislobased_Nslip(matID) + +!* Dislocation density evolution + constitutive_dislobased_dotState = 0.0_pReal + do i = 1,n + tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) + 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)*& + 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))/& + (constitutive_dislobased_c4(matID)*constitutive_dislobased_bg(matID)) + + 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))/& + ! (kB*Temperature))*exp(-constitutive_dislobased_Qsd(matID)/(kB*Temperature)) + + constitutive_dislobased_dotState(i) = locks - athermal_recovery + endif + enddo + return end function +!**************************************************************** +!* calculates the rate of change of temperature * +!**************************************************************** +pure function constitutive_dislobased_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el) + + !*** variables and functions from other modules ***! + use prec, only: pReal,pInt,p_vec + use lattice, only: lattice_Sslip_v + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance + implicit none + + !*** input variables ***! + real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: Temperature + integer(pInt), intent(in):: ipc, & ! grain number + ip, & ! integration point number + el ! element number + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure + + !*** output variables ***! + real(pReal) constitutive_dislobased_dotTemperature ! rate of change of temparature + + ! calculate dotTemperature + constitutive_dislobased_dotTemperature = 0.0_pReal + + return +endfunction + + + pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) !********************************************************************* !* return array of constitutive results * @@ -729,9 +587,10 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i !* - ip : current integration point * !* - el : current element * !********************************************************************* - use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_shearTwin - use mesh, only: mesh_NcpElems,mesh_maxNips + use prec, only: pReal,pInt,p_vec + use math, only: math_mul6x6 + use lattice, only: lattice_Sslip_v + use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput implicit none @@ -740,85 +599,37 @@ pure function constitutive_dislobased_postResults(Tstar_v,Temperature,dt,state,i real(pReal), intent(in) :: dt,Temperature real(pReal), dimension(6), intent(in) :: Tstar_v type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state - integer(pInt) matID,o,i,c,ns,nt - real(pReal) sumf,tau_slip,tau_twin + integer(pInt) matID,o,i,c,n + real(pReal) tau_slip, active_rate real(pReal), dimension(constitutive_dislobased_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & constitutive_dislobased_postResults - !* Shortened notation matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - ns = constitutive_dislobased_Nslip(matID) - nt = constitutive_dislobased_Ntwin(matID) - - !* Total twin volume fraction - sumf = 0.0_pReal - if (nt > 0_pInt) sumf = sum(state(ipc,ip,el)%p((ns+1):(ns+nt))) - - !* Required output + n = constitutive_dislobased_Nslip(matID) c = 0_pInt constitutive_dislobased_postResults = 0.0_pReal do o = 1,phase_Noutput(material_phase(ipc,ip,el)) select case(constitutive_dislobased_output(o,matID)) - - case ('state_slip') - constitutive_dislobased_postResults(c+1:c+ns) = state(ipc,ip,el)%p(1:ns) - c = c + ns - - case ('rateofshear_slip') - do i = 1,ns - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) - - if ( abs(tau_slip) > state(ipc,ip,el)%p(6*ns+3*nt+i) ) then - constitutive_dislobased_postResults(c+i) = state(ipc,ip,el)%p(9*ns+5*nt+i)*sign(1.0_pReal,tau_slip)*& - sinh(((abs(tau_slip)-state(ipc,ip,el)%p(6*ns+3*nt+i))*state(ipc,ip,el)%p(7*ns+4*nt+i))/(kB*Temperature) ) + case ('dislodensity') + constitutive_dislobased_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n) + c = c + n + case ('rateofshear') + do i = 1,n + tau_slip = math_mul6x6(Tstar_v,lattice_Sslip_v(:,i,constitutive_dislobased_structure(matID))) + if ((abs(tau_slip)-state(ipc,ip,el)%p(3*n+i))>0) then + 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 + constitutive_dislobased_postResults(c+i) = 0.0_pReal + endif enddo - c = c + ns - - case ('mfp_slip') - constitutive_dislobased_postResults(c+1:c+ns) = state(ipc,ip,el)%p((5*ns+2*nt+1):(6*ns+2*nt)) - c = c + ns - - case ('thresholdstress_slip') - constitutive_dislobased_postResults(c+1:c+ns) = state(ipc,ip,el)%p((6*ns+3*nt+1):(7*ns+3*nt)) - c = c + ns - - case ('state_twin') - if (nt > 0_pInt) constitutive_dislobased_postResults(c+1:c+nt) = state(ipc,ip,el)%p((ns+1):(ns+nt)) - c = c + nt - - case ('rateofshear_twin') - if (nt > 0_pInt) then - do i = 1,nt - tau_twin = dot_product(Tstar_v,lattice_Stwin_v(:,i,constitutive_dislobased_structure(matID))) - - if ( tau_twin > 0.0_pReal ) then - constitutive_dislobased_postResults(c+i) = (constitutive_dislobased_fmax(matID) - sumf)*& - lattice_shearTwin(i,constitutive_dislobased_structure(matID))*& - state(ipc,ip,el)%p(8*ns+4*nt+i)*constitutive_dislobased_Ndot0(matID)*& - exp(-(state(ipc,ip,el)%p(7*ns+3*nt+i)/tau_twin)**10.0_pReal) - else - constitutive_dislobased_postResults(c+i) = 0.0_pReal - endif - enddo - endif - c = c + nt - - case ('mfp_twin') - if (nt > 0_pInt) constitutive_dislobased_postResults(c+1:c+nt) = state(ipc,ip,el)%p((6*ns+2*nt+1):(6*ns+3*nt)) - c = c + nt - - case ('thresholdstress_twin') - if (nt > 0_pInt) constitutive_dislobased_postResults(c+1:c+nt) = state(ipc,ip,el)%p((7*ns+3*nt+1):(7*ns+4*nt)) - c = c + nt - + c = c + n end select enddo return + end function END MODULE \ No newline at end of file diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 06550c7e2..5d6360544 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -18,7 +18,7 @@ ! gdot0 0.001 ! n 20 ! h0 75e6 -! tau_sat 63e6 +! tausat 63e6 ! w0 2.25 MODULE constitutive_j2 @@ -32,7 +32,8 @@ MODULE constitutive_j2 integer(pInt), dimension(:), allocatable :: constitutive_j2_sizeDotState, & constitutive_j2_sizeState, & constitutive_j2_sizePostResults - character(len=64), dimension(:,:), allocatable :: constitutive_j2_output + integer(pInt), dimension(:,:), allocatable,target :: constitutive_j2_sizePostResult ! size of each post result output + character(len=64), dimension(:,:), allocatable,target :: constitutive_j2_output ! name of each post result output real(pReal), dimension(:), allocatable :: constitutive_j2_C11 real(pReal), dimension(:), allocatable :: constitutive_j2_C12 real(pReal), dimension(:,:,:), allocatable :: constitutive_j2_Cslip_66 @@ -42,7 +43,7 @@ MODULE constitutive_j2 real(pReal), dimension(:), allocatable :: constitutive_j2_gdot0 real(pReal), dimension(:), allocatable :: constitutive_j2_n real(pReal), dimension(:), allocatable :: constitutive_j2_h0 - real(pReal), dimension(:), allocatable :: constitutive_j2_tau_sat + real(pReal), dimension(:), allocatable :: constitutive_j2_tausat real(pReal), dimension(:), allocatable :: constitutive_j2_w0 @@ -69,9 +70,13 @@ subroutine constitutive_j2_init(file) integer(pInt), intent(in) :: file integer(pInt), parameter :: maxNchunks = 7 integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) section, maxNinstance, i,j,k,l, output + integer(pInt) section, maxNinstance, i,j,k,l, output, mySize character(len=64) tag character(len=1024) line + + write(6,*) + write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_j2_label,' init -+>>>' + write(6,*) maxNinstance = count(phase_constitution == constitutive_j2_label) if (maxNinstance == 0) return @@ -79,6 +84,8 @@ subroutine constitutive_j2_init(file) allocate(constitutive_j2_sizeDotState(maxNinstance)) ; constitutive_j2_sizeDotState = 0_pInt allocate(constitutive_j2_sizeState(maxNinstance)) ; constitutive_j2_sizeState = 0_pInt allocate(constitutive_j2_sizePostResults(maxNinstance)); constitutive_j2_sizePostResults = 0_pInt + allocate(constitutive_j2_sizePostResult(maxval(phase_Noutput), & + maxNinstance)) ; constitutive_j2_sizePostResult = 0_pInt allocate(constitutive_j2_output(maxval(phase_Noutput), & maxNinstance)) ; constitutive_j2_output = '' allocate(constitutive_j2_C11(maxNinstance)) ; constitutive_j2_C11 = 0.0_pReal @@ -89,7 +96,7 @@ subroutine constitutive_j2_init(file) allocate(constitutive_j2_gdot0(maxNinstance)) ; constitutive_j2_gdot0 = 0.0_pReal allocate(constitutive_j2_n(maxNinstance)) ; constitutive_j2_n = 0.0_pReal allocate(constitutive_j2_h0(maxNinstance)) ; constitutive_j2_h0 = 0.0_pReal - allocate(constitutive_j2_tau_sat(maxNinstance)) ; constitutive_j2_tau_sat = 0.0_pReal + allocate(constitutive_j2_tausat(maxNinstance)) ; constitutive_j2_tausat = 0.0_pReal allocate(constitutive_j2_w0(maxNinstance)) ; constitutive_j2_w0 = 0.0_pReal rewind(file) @@ -128,8 +135,8 @@ subroutine constitutive_j2_init(file) constitutive_j2_n(i) = IO_floatValue(line,positions,2) case ('h0') constitutive_j2_h0(i) = IO_floatValue(line,positions,2) - case ('tau_sat') - constitutive_j2_tau_sat(i) = IO_floatValue(line,positions,2) + case ('tausat') + constitutive_j2_tausat(i) = IO_floatValue(line,positions,2) case ('w0') constitutive_j2_w0(i) = IO_floatValue(line,positions,2) case ('taylorfactor') @@ -139,35 +146,40 @@ subroutine constitutive_j2_init(file) enddo 100 do i = 1,maxNinstance ! sanity checks - if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(203) - if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(204) - if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(205) - if (constitutive_j2_h0(i) <= 0.0_pReal) call IO_error(206) - if (constitutive_j2_tau_sat(i) <= 0.0_pReal) call IO_error(207) - if (constitutive_j2_w0(i) <= 0.0_pReal) call IO_error(208) + if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(210) + if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211) + if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(212) + if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(213) + if (constitutive_j2_w0(i) <= 0.0_pReal) call IO_error(241) if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240) enddo do i = 1,maxNinstance + do j = 1,maxval(phase_Noutput) + select case(constitutive_j2_output(j,i)) + case('flowstress') + mySize = 1_pInt + case('strainrate') + mySize = 1_pInt + case default + mySize = 0_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + constitutive_j2_sizePostResult(j,i) = mySize + constitutive_j2_sizePostResults(i) = & + constitutive_j2_sizePostResults(i) + mySize + endif + enddo + constitutive_j2_sizeDotState(i) = 1 constitutive_j2_sizeState(i) = 1 - do j = 1,maxval(phase_Noutput) - select case(constitutive_j2_output(j,i)) - case('flowstress') - constitutive_j2_sizePostResults(i) = & - constitutive_j2_sizePostResults(i) + 1 - case('strainrate') - constitutive_j2_sizePostResults(i) = & - constitutive_j2_sizePostResults(i) + 1 - end select - enddo - forall(k=1:3) forall(j=1:3) & - constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i) - constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i) - constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i)) + constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i) + constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i) + constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i)) end forall constitutive_j2_Cslip_66(:,:,i) = & math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(:,:,i))) @@ -184,23 +196,12 @@ endsubroutine !********************************************************************* pure function constitutive_j2_stateInit(myInstance) - !*** variables and functions from other modules ***! use prec, only: pReal,pInt - implicit none - !*** input variables ***! integer(pInt), intent(in) :: myInstance - - !*** output variables ***! real(pReal), dimension(1) :: constitutive_j2_stateInit - !*** local variables ***! - - !*** global variables ***! - ! constitutive_j2_tau0 - - constitutive_j2_stateInit = constitutive_j2_tau0(myInstance) return @@ -221,7 +222,6 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el) use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance implicit none -!* Definition of variables integer(pInt), intent(in) :: ipc,ip,el integer(pInt) matID real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC @@ -305,12 +305,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, l, & m, & n - - !*** global variables ***! - ! constitutive_j2_gdot0 - ! constitutive_j2_fTaylor - ! constitutive_j2_n - matID = phase_constitutionInstance(material_phase(g,ip,el)) @@ -384,66 +378,58 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el) norm_Tstar_dev ! euclidean norm of Tstar_dev integer(pInt) matID - !*** global variables ***! - ! constitutive_j2_gdot0 - ! constitutive_j2_fTaylor - ! constitutive_j2_n - ! constitutive_j2_h0 - ! constitutive_j2_tau_sat - ! constitutive_j2_w0 - matID = phase_constitutionInstance(material_phase(g,ip,el)) - ! caclulate deviatoric part of 2nd Piola-Kirchhoff stress - Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3_pReal + ! deviatoric part of 2nd Piola-Kirchhoff stress + Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal Tstar_dev_v(4:6) = Tstar_v(4:6) norm_Tstar_dev = dsqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) - ! Calculation of gamma_dot + ! gamma_dot gamma_dot = constitutive_j2_gdot0(matID) * ( dsqrt(1.5_pReal) * norm_Tstar_dev & / &!--------------------------------------------------- (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID) - ! calculate hardening coefficient + ! hardening coefficient hardening = constitutive_j2_h0(matID) * & - ( 1.0_pReal - state(g,ip,el)%p(1) / constitutive_j2_tau_sat(matID) ) ** constitutive_j2_w0(matID) + ( 1.0_pReal - state(g,ip,el)%p(1) / constitutive_j2_tausat(matID) ) ** constitutive_j2_w0(matID) - ! finally calculate dotState + ! dotState constitutive_j2_dotState = hardening * gamma_dot return -endfunction - - -!**************************************************************** -!* calculates the rate of change of microstructure * -!**************************************************************** -pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el) - - !*** variables and functions from other modules ***! - use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance - implicit none - - !*** input variables ***! - real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: Temperature - integer(pInt), intent(in):: g, & ! grain number - ip, & ! integration point number - el ! element number - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure - - !*** output variables ***! - real(pReal) constitutive_j2_dotTemperature ! evolution of state variable - - ! calculate dotTemperature - constitutive_j2_dotTemperature = 0.0_pReal - - return +endfunction + + +!**************************************************************** +!* calculates the rate of change of temperature * +!**************************************************************** +pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el) + + !*** variables and functions from other modules ***! + use prec, only: pReal,pInt,p_vec + use lattice, only: lattice_Sslip_v + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance + implicit none + + !*** input variables ***! + real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: Temperature + integer(pInt), intent(in):: g, & ! grain number + ip, & ! integration point number + el ! element number + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure + + !*** output variables ***! + real(pReal) constitutive_j2_dotTemperature ! rate of change of temperature + + ! calculate dotTemperature + constitutive_j2_dotTemperature = 0.0_pReal + + return endfunction diff --git a/code/constitutive_phenomenological.f90 b/code/constitutive_phenomenological.f90 deleted file mode 100644 index 7fe2f6581..000000000 --- a/code/constitutive_phenomenological.f90 +++ /dev/null @@ -1,519 +0,0 @@ - -!***************************************************** -!* Module: CONSTITUTIVE_PHENOMENOLOGICAL * -!***************************************************** -!* contains: * -!* - constitutive equations * -!* - parameters definition * -!***************************************************** - -! [Alu] -! constitution phenomenological -! (output) slipresistance -! (output) rateofshear -! lattice_structure 1 -! Nslip 12 -! -! c11 106.75e9 -! c12 60.41e9 -! c44 28.34e9 -! -! s0_slip 31e6 -! gdot0_slip 0.001 -! n_slip 20 -! h0 75e6 -! s_sat 63e6 -! w0 2.25 -! latent_ratio 1.4 - -MODULE constitutive_phenomenological - -!*** Include other modules *** - use prec, only: pReal,pInt - implicit none - - character (len=*), parameter :: constitutive_phenomenological_label = 'phenomenological' - - integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_sizeDotState, & - constitutive_phenomenological_sizeState, & - constitutive_phenomenological_sizePostResults - 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_Nslip - - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_CoverA - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C11 - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C12 - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C13 - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C33 - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C44 - real(pReal), dimension(:,:,:), allocatable :: constitutive_phenomenological_Cslip_66 -!* Visco-plastic constitutive_phenomenological parameters - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_s0_slip - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_gdot0_slip - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_n_slip - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_h0 - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_s_sat - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_w0 - real(pReal), dimension(:), allocatable :: constitutive_phenomenological_latent - real(pReal), dimension(:,:,:), allocatable :: constitutive_phenomenological_HardeningMatrix - - - -CONTAINS -!**************************************** -!* - constitutive_init -!* - constitutive_stateInit -!* - constitutive_homogenizedC -!* - constitutive_microstructure -!* - constitutive_LpAndItsTangent -!* - consistutive_dotState -!* - consistutive_postResults -!**************************************** - - -subroutine constitutive_phenomenological_init(file) -!************************************** -!* Module initialization * -!************************************** - use prec, only: pInt, pReal - use math, only: math_Mandel3333to66, math_Voigt66to3333 - use IO - use material - - use lattice, only: lattice_initializeStructure - integer(pInt), intent(in) :: file - integer(pInt), parameter :: maxNchunks = 7 - integer(pInt), dimension(1+2*maxNchunks) :: positions - integer(pInt) section, maxNinstance, i,j,k, output - character(len=64) tag - character(len=1024) line - - maxNinstance = count(phase_constitution == constitutive_phenomenological_label) - if (maxNinstance == 0) return - - allocate(constitutive_phenomenological_sizeDotState(maxNinstance)) ; constitutive_phenomenological_sizeDotState = 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_output(maxval(phase_Noutput), & - 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_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_C12(maxNinstance)) ; constitutive_phenomenological_C12 = 0.0_pReal - allocate(constitutive_phenomenological_C13(maxNinstance)) ; constitutive_phenomenological_C13 = 0.0_pReal - allocate(constitutive_phenomenological_C33(maxNinstance)) ; constitutive_phenomenological_C33 = 0.0_pReal - allocate(constitutive_phenomenological_C44(maxNinstance)) ; constitutive_phenomenological_C44 = 0.0_pReal - allocate(constitutive_phenomenological_Cslip_66(6,6,maxNinstance)) ; constitutive_phenomenological_Cslip_66 = 0.0_pReal - allocate(constitutive_phenomenological_s0_slip(maxNinstance)) ; constitutive_phenomenological_s0_slip = 0.0_pReal - allocate(constitutive_phenomenological_gdot0_slip(maxNinstance)) ; constitutive_phenomenological_gdot0_slip = 0.0_pReal - allocate(constitutive_phenomenological_n_slip(maxNinstance)) ; constitutive_phenomenological_n_slip = 0.0_pReal - allocate(constitutive_phenomenological_h0(maxNinstance)) ; constitutive_phenomenological_h0 = 0.0_pReal - allocate(constitutive_phenomenological_s_sat(maxNinstance)) ; constitutive_phenomenological_s_sat = 0.0_pReal - allocate(constitutive_phenomenological_w0(maxNinstance)) ; constitutive_phenomenological_w0 = 0.0_pReal - allocate(constitutive_phenomenological_latent(maxNinstance)) ; constitutive_phenomenological_latent = 0.0_pReal - - rewind(file) - line = '' - section = 0 - - do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to - read(file,'(a1024)',END=100) line - enddo - - do ! read thru sections of phase part - read(file,'(a1024)',END=100) line - if (IO_isBlank(line)) cycle ! skip empty lines - if (IO_getTag(line,'<','>') /= '') exit ! stop at next part - if (IO_getTag(line,'[',']') /= '') then ! next section - section = section + 1 - output = 0 ! reset output counter - endif - if (section > 0 .and. phase_constitution(section) == constitutive_phenomenological_label) then ! one of my sections - i = phase_constitutionInstance(section) ! which instance of my constitution is present phase - positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key - select case(tag) - case ('(output)') - output = output + 1 - constitutive_phenomenological_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) - case ('lattice_structure') - constitutive_phenomenological_structureName(i) = IO_lc(IO_stringValue(line,positions,2)) - case ('nslip') - constitutive_phenomenological_Nslip(i) = IO_intValue(line,positions,2) - case ('covera_ratio') - constitutive_phenomenological_CoverA(i) = IO_floatValue(line,positions,2) - case ('c11') - constitutive_phenomenological_C11(i) = IO_floatValue(line,positions,2) - case ('c12') - constitutive_phenomenological_C12(i) = IO_floatValue(line,positions,2) - case ('c13') - constitutive_phenomenological_C13(i) = IO_floatValue(line,positions,2) - case ('c33') - constitutive_phenomenological_C33(i) = IO_floatValue(line,positions,2) - case ('c44') - constitutive_phenomenological_C44(i) = IO_floatValue(line,positions,2) - case ('s0_slip') - constitutive_phenomenological_s0_slip(i) = IO_floatValue(line,positions,2) - case ('gdot0_slip') - constitutive_phenomenological_gdot0_slip(i) = IO_floatValue(line,positions,2) - case ('n_slip') - constitutive_phenomenological_n_slip(i) = IO_floatValue(line,positions,2) - case ('h0') - constitutive_phenomenological_h0(i) = IO_floatValue(line,positions,2) - case ('s_sat') - constitutive_phenomenological_s_sat(i) = IO_floatValue(line,positions,2) - case ('w0') - constitutive_phenomenological_w0(i) = IO_floatValue(line,positions,2) - case ('latent_ratio') - constitutive_phenomenological_latent(i) = IO_floatValue(line,positions,2) - end select - endif - enddo - -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. & - constitutive_phenomenological_structure(i) > 3) call IO_error(201) - if (constitutive_phenomenological_Nslip(i) < 1) call IO_error(202) - if (constitutive_phenomenological_s0_slip(i) < 0.0_pReal) call IO_error(203) - if (constitutive_phenomenological_gdot0_slip(i) <= 0.0_pReal) call IO_error(204) - if (constitutive_phenomenological_n_slip(i) <= 0.0_pReal) call IO_error(205) - if (constitutive_phenomenological_h0(i) <= 0.0_pReal) call IO_error(206) - if (constitutive_phenomenological_s_sat(i) <= 0.0_pReal) call IO_error(207) - if (constitutive_phenomenological_w0(i) <= 0.0_pReal) call IO_error(208) - if (constitutive_phenomenological_latent(i) <= 0.0_pReal) call IO_error(209) - enddo - - allocate(constitutive_phenomenological_hardeningMatrix(maxval(constitutive_phenomenological_Nslip),& - maxval(constitutive_phenomenological_Nslip),& - maxNinstance)) - - do i = 1,maxNinstance - constitutive_phenomenological_sizeDotState(i) = constitutive_phenomenological_Nslip(i) - constitutive_phenomenological_sizeState(i) = constitutive_phenomenological_Nslip(i) - - do j = 1,maxval(phase_Noutput) - select case(constitutive_phenomenological_output(j,i)) - case('slipresistance') - constitutive_phenomenological_sizePostResults(i) = & - constitutive_phenomenological_sizePostResults(i) + constitutive_phenomenological_Nslip(i) - case('rateofshear') - constitutive_phenomenological_sizePostResults(i) = & - constitutive_phenomenological_sizePostResults(i) + constitutive_phenomenological_Nslip(i) - end select - enddo - - select case (constitutive_phenomenological_structure(i)) - case(1:2) ! cubic(s) - forall(k=1:3) - forall(j=1:3) & - constitutive_phenomenological_Cslip_66(k,j,i) = constitutive_phenomenological_C12(i) - constitutive_phenomenological_Cslip_66(k,k,i) = constitutive_phenomenological_C11(i) - constitutive_phenomenological_Cslip_66(k+3,k+3,i) = constitutive_phenomenological_C44(i) - end forall - case(3) ! hcp - constitutive_phenomenological_Cslip_66(1,1,i) = constitutive_phenomenological_C11(i) - constitutive_phenomenological_Cslip_66(2,2,i) = constitutive_phenomenological_C11(i) - constitutive_phenomenological_Cslip_66(3,3,i) = constitutive_phenomenological_C33(i) - constitutive_phenomenological_Cslip_66(1,2,i) = constitutive_phenomenological_C12(i) - constitutive_phenomenological_Cslip_66(2,1,i) = constitutive_phenomenological_C12(i) - constitutive_phenomenological_Cslip_66(1,3,i) = constitutive_phenomenological_C13(i) - constitutive_phenomenological_Cslip_66(3,1,i) = constitutive_phenomenological_C13(i) - constitutive_phenomenological_Cslip_66(2,3,i) = constitutive_phenomenological_C13(i) - constitutive_phenomenological_Cslip_66(3,2,i) = constitutive_phenomenological_C13(i) - constitutive_phenomenological_Cslip_66(4,4,i) = constitutive_phenomenological_C44(i) - constitutive_phenomenological_Cslip_66(5,5,i) = constitutive_phenomenological_C44(i) - constitutive_phenomenological_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_phenomenological_C11(i)- & - constitutive_phenomenological_C12(i)) - end select - constitutive_phenomenological_Cslip_66(:,:,i) = & - math_Mandel3333to66(math_Voigt66to3333(constitutive_phenomenological_Cslip_66(:,:,i))) - - constitutive_phenomenological_hardeningMatrix(:,:,i) = constitutive_phenomenological_latent(i) - forall (j = 1:constitutive_phenomenological_Nslip(i)) & - constitutive_phenomenological_hardeningMatrix(j,j,i) = 1.0_pReal - enddo - - return - -end subroutine - - -function constitutive_phenomenological_stateInit(myInstance) -!********************************************************************* -!* initial microstructural state * -!********************************************************************* - use prec, only: pReal,pInt - implicit none - -!* Definition of variables - integer(pInt), intent(in) :: myInstance - real(pReal), dimension(constitutive_phenomenological_Nslip(myInstance)) :: constitutive_phenomenological_stateInit - - constitutive_phenomenological_stateInit = constitutive_phenomenological_s0_slip(myInstance) - - return -end function - - -function constitutive_phenomenological_homogenizedC(state,ipc,ip,el) -!********************************************************************* -!* homogenized elacticity matrix * -!* INPUT: * -!* - state : state variables * -!* - ipc : component-ID of current integration point * -!* - ip : current integration point * -!* - el : current element * -!********************************************************************* - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance - implicit none - -!* Definition of variables - integer(pInt), intent(in) :: ipc,ip,el - integer(pInt) matID - real(pReal), dimension(6,6) :: constitutive_phenomenological_homogenizedC - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - - matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - constitutive_phenomenological_homogenizedC = constitutive_phenomenological_Cslip_66(:,:,matID) - - return - -end function - - -subroutine constitutive_phenomenological_microstructure(Temperature,state,ipc,ip,el) -!********************************************************************* -!* calculate derived quantities from state (not used here) * -!* INPUT: * -!* - Tp : temperature * -!* - ipc : component-ID of current integration point * -!* - ip : current integration point * -!* - el : current element * -!********************************************************************* - use prec, only: pReal,pInt,p_vec - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance - implicit none - -!* Definition of variables - integer(pInt) ipc,ip,el, matID - real(pReal) Temperature - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - - matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - -end subroutine - - -subroutine constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el) -!********************************************************************* -!* plastic velocity gradient and its tangent * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - ipc : component-ID at current integration point * -!* - ip : current integration point * -!* - el : current element * -!* OUTPUT: * -!* - Lp : plastic velocity gradient * -!* - dLp_dTstar : derivative of Lp (4th-rank tensor) * -!********************************************************************* - use prec, only: pReal,pInt,p_vec - use math, only: math_Plain3333to99 - use lattice, only: lattice_Sslip,lattice_Sslip_v - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance - - implicit none - -!* Definition of variables - integer(pInt) ipc,ip,el - integer(pInt) matID,i,k,l,m,n - real(pReal) Temperature - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - real(pReal), dimension(6) :: Tstar_v - real(pReal), dimension(3,3) :: Lp - real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 - real(pReal), dimension(9,9) :: dLp_dTstar - real(pReal), dimension(constitutive_phenomenological_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & - gdot_slip,dgdot_dtauslip,tau_slip - - matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - -!* Calculation of Lp - Lp = 0.0_pReal - do i = 1,constitutive_phenomenological_Nslip(matID) - tau_slip(i) = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) - - gdot_slip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**& - constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip(i)) - - Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_phenomenological_structure(matID)) - enddo - -!* Calculation of the tangent of Lp - dLp_dTstar3333 = 0.0_pReal - dLp_dTstar = 0.0_pReal - do i = 1,constitutive_phenomenological_Nslip(matID) - - dgdot_dtauslip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**& - (constitutive_phenomenological_n_slip(matID)-1.0_pReal)*& - constitutive_phenomenological_n_slip(matID)/state(ipc,ip,el)%p(i) - - 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) + & - dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_phenomenological_structure(matID))* & - lattice_Sslip(m,n,i,constitutive_phenomenological_structure(matID)) - enddo - dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) - - return -end subroutine - - -function constitutive_phenomenological_dotState(Tstar_v,Temperature,state,ipc,ip,el) -!********************************************************************* -!* rate of change of microstructure * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - ipc : component-ID at current integration point * -!* - ip : current integration point * -!* - el : current element * -!* OUTPUT: * -!* - constitutive_dotState : evolution of state variable * -!********************************************************************* - use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance - implicit none - -!* Definition of variables - integer(pInt) ipc,ip,el - integer(pInt) matID,i,n - real(pReal) Temperature,tau_slip,gdot_slip - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - real(pReal), dimension(6) :: Tstar_v - real(pReal), dimension(constitutive_phenomenological_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & - constitutive_phenomenological_dotState,self_hardening - - matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - n = constitutive_phenomenological_Nslip(matID) - -!* Self-Hardening of each system - do i = 1,n - - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) - - gdot_slip = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip)/state(ipc,ip,el)%p(i))**& - constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip) - - self_hardening(i) = constitutive_phenomenological_h0(matID)*(1.0_pReal-state(ipc,ip,el)%p(i)/& - constitutive_phenomenological_s_sat(matID))**constitutive_phenomenological_w0(matID)*abs(gdot_slip) - enddo - -!$OMP CRITICAL (evilmatmul) - constitutive_phenomenological_dotState = matmul(constitutive_phenomenological_hardeningMatrix(1:n,1:n,matID),self_hardening) -!$OMP END CRITICAL (evilmatmul) - return - -end function - - -function constitutive_phenomenological_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el) -!********************************************************************* -!* rate of change of microstructure * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - ipc : component-ID at current integration point * -!* - ip : current integration point * -!* - el : current element * -!* OUTPUT: * -!* - constitutive_dotTemperature : evolution of temperature * -!********************************************************************* - use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance - implicit none - -!* Definition of variables - integer(pInt) ipc,ip,el - integer(pInt) matID,i,n - real(pReal) Temperature - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state - real(pReal), dimension(6) :: Tstar_v - real(pReal) constitutive_phenomenological_dotTemperature - - constitutive_phenomenological_dotTemperature = 0.0_pReal - - return -end function - - -pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) -!********************************************************************* -!* return array of constitutive results * -!* INPUT: * -!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * -!* - dt : current time increment * -!* - ipc : component-ID at current integration point * -!* - ip : current integration point * -!* - el : current element * -!********************************************************************* - use prec, only: pReal,pInt,p_vec - use lattice, only: lattice_Sslip_v - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput - implicit none - -!* Definition of variables - integer(pInt), intent(in) :: ipc,ip,el - real(pReal), intent(in) :: dt,Temperature - real(pReal), dimension(6), intent(in) :: Tstar_v - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state - integer(pInt) matID,o,i,c,n - real(pReal) tau_slip - real(pReal), dimension(constitutive_phenomenological_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & - constitutive_phenomenological_postResults - - matID = phase_constitutionInstance(material_phase(ipc,ip,el)) - n = constitutive_phenomenological_Nslip(matID) - c = 0_pInt - constitutive_phenomenological_postResults = 0.0_pReal - - do o = 1,phase_Noutput(material_phase(ipc,ip,el)) - select case(constitutive_phenomenological_output(o,matID)) - - case ('slipresistance') - constitutive_phenomenological_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n) - c = c + n - - case ('rateofshear') - do i = 1,n - tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID))) - constitutive_phenomenological_postResults(c+i) = sign(1.0_pReal,tau_slip)*constitutive_phenomenological_gdot0_slip(matID)*& - (abs(tau_slip)/state(ipc,ip,el)%p(i))**& - constitutive_phenomenological_n_slip(matID) - enddo - c = c + n - - end select - enddo - - return - -end function - -END MODULE diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 new file mode 100644 index 000000000..45646d0af --- /dev/null +++ b/code/constitutive_phenopowerlaw.f90 @@ -0,0 +1,912 @@ + +!***************************************************** +!* Module: CONSTITUTIVE_PHENOPOWERLAW * +!***************************************************** +!* contains: * +!* - constitutive equations * +!* - parameters definition * +!***************************************************** + +! [Alu] +! constitution phenopowerlaw +! (output) resistance_slip +! (output) shearrate_slip +! (output) resolvedstress_slip +! (output) totalshear +! (output) resistance_twin +! (output) shearrate_twin +! (output) resolvedstress_twin +! (output) totalvolfrac +! lattice_structure hex +! covera_ratio 1.57 +! Nslip 3 3 6 12 # per family +! Ntwin 6 6 6 6 # per family +! +! c11 106.75e9 +! c12 60.41e9 +! c44 28.34e9 +! +! gdot0_slip 0.001 +! n_slip 20 +! tau0_slip 31e6 31e6 60e6 123e6 # per family +! tausat_slip 63e6 90e6 200e6 400e6 # per family +! gdot0_twin 0.001 +! n_twin 20 +! tau0_twin 31e6 31e6 60e6 123e6 # per family +! s_pr 100e6 # push-up factor for slip saturation due to twinning +! twin_b 2 +! twin_c 25 +! twin_d 6 +! twin_e 9 +! h0_slipslip 75e6 +! h0_sliptwin 75e6 +! h0_twinslip 75e6 +! h0_twintwin 75e6 +! interaction_slipslip 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 +! interaction_sliptwin 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 +! interaction_twinslip 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 +! interaction_twintwin 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 + +MODULE constitutive_phenopowerlaw + +!*** Include other modules *** + use prec, only: pReal,pInt + implicit none + + character (len=*), parameter :: constitutive_phenopowerlaw_label = 'phenopowerlaw' + + integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_sizeDotState, & + constitutive_phenopowerlaw_sizeState, & + constitutive_phenopowerlaw_sizePostResults ! cumulative size of post results + integer(pInt), dimension(:,:), allocatable,target :: constitutive_phenopowerlaw_sizePostResult ! size of each post result output + character(len=64), dimension(:,:), allocatable,target :: constitutive_phenopowerlaw_output ! name of each post result output + + character(len=32), dimension(:), allocatable :: constitutive_phenopowerlaw_structureName + integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_structure + integer(pInt), dimension(:,:), allocatable :: constitutive_phenopowerlaw_Nslip ! active number of slip systems per family + integer(pInt), dimension(:,:), allocatable :: constitutive_phenopowerlaw_Ntwin ! active number of twin systems per family + integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_totalNslip ! no. of slip system used in simulation + integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_totalNtwin ! no. of twin system used in simulation + + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_CoverA + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C11 + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C12 + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C13 + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C33 + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C44 + real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_Cslip_66 +!* Visco-plastic constitutive_phenomenological parameters + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_gdot0_slip + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_n_slip + real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tau0_slip + real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tausat_slip + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_gdot0_twin + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_n_twin + real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tau0_twin + + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_spr + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinB + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinC + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinD + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinE + + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_slipslip + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_sliptwin + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_twinslip + real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_twintwin + + real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_slipslip + real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_sliptwin + real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_twinslip + real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_twintwin + + real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_slipslip + real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_sliptwin + real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twinslip + real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twintwin + + +CONTAINS +!**************************************** +!* - constitutive_init +!* - constitutive_stateInit +!* - constitutive_homogenizedC +!* - constitutive_microstructure +!* - constitutive_LpAndItsTangent +!* - consistutive_dotState +!* - consistutive_postResults +!**************************************** + + +subroutine constitutive_phenopowerlaw_init(file) +!************************************** +!* Module initialization * +!************************************** + use prec, only: pInt, pReal + use math, only: math_Mandel3333to66, math_Voigt66to3333 + use IO + use material + + use lattice, only: lattice_initializeStructure, lattice_maxNslipFamily, lattice_maxNtwinFamily, & + lattice_maxNinteraction, lattice_NslipSystem, lattice_NtwinSystem, & + lattice_interactionSlipSlip,lattice_interactionSlipTwin,lattice_interactionTwinSlip,lattice_interactionTwinTwin + integer(pInt), intent(in) :: file + integer(pInt), parameter :: maxNchunks = 21 + integer(pInt), dimension(1+2*maxNchunks) :: positions + integer(pInt) section, maxNinstance, i,j,k,l,m, output, mySize + character(len=64) tag + character(len=1024) line + + write(6,*) + write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_phenopowerlaw_label,' init -+>>>' + write(6,*) + + + maxNinstance = count(phase_constitution == constitutive_phenopowerlaw_label) + if (maxNinstance == 0) return + + allocate(constitutive_phenopowerlaw_sizeDotState(maxNinstance)) ; constitutive_phenopowerlaw_sizeDotState = 0_pInt + allocate(constitutive_phenopowerlaw_sizeState(maxNinstance)) ; constitutive_phenopowerlaw_sizeState = 0_pInt + allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance)); constitutive_phenopowerlaw_sizePostResults = 0_pInt + allocate(constitutive_phenopowerlaw_sizePostResult(maxval(phase_Noutput), & + maxNinstance)) ; constitutive_phenopowerlaw_sizePostResult = 0_pInt + allocate(constitutive_phenopowerlaw_output(maxval(phase_Noutput), & + maxNinstance)) ; constitutive_phenopowerlaw_output = '' + + allocate(constitutive_phenopowerlaw_structureName(maxNinstance)) ; constitutive_phenopowerlaw_structureName = '' + allocate(constitutive_phenopowerlaw_structure(maxNinstance)) ; constitutive_phenopowerlaw_structure = 0_pInt + allocate(constitutive_phenopowerlaw_Nslip(lattice_maxNslipFamily,& + maxNinstance)) ; constitutive_phenopowerlaw_Nslip = 0_pInt + allocate(constitutive_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,& + maxNinstance)) ; constitutive_phenopowerlaw_Ntwin = 0_pInt + + allocate(constitutive_phenopowerlaw_totalNslip(maxNinstance)) ; constitutive_phenopowerlaw_totalNslip = 0_pInt !no. of slip system used in simulation (YJ.RO) + allocate(constitutive_phenopowerlaw_totalNtwin(maxNinstance)) ; constitutive_phenopowerlaw_totalNtwin = 0_pInt !no. of twin system used in simulation (YJ.RO) + + allocate(constitutive_phenopowerlaw_CoverA(maxNinstance)) ; constitutive_phenopowerlaw_CoverA = 0.0_pReal + allocate(constitutive_phenopowerlaw_C11(maxNinstance)) ; constitutive_phenopowerlaw_C11 = 0.0_pReal + allocate(constitutive_phenopowerlaw_C12(maxNinstance)) ; constitutive_phenopowerlaw_C12 = 0.0_pReal + allocate(constitutive_phenopowerlaw_C13(maxNinstance)) ; constitutive_phenopowerlaw_C13 = 0.0_pReal + allocate(constitutive_phenopowerlaw_C33(maxNinstance)) ; constitutive_phenopowerlaw_C33 = 0.0_pReal + allocate(constitutive_phenopowerlaw_C44(maxNinstance)) ; constitutive_phenopowerlaw_C44 = 0.0_pReal + allocate(constitutive_phenopowerlaw_Cslip_66(6,6,maxNinstance)) ; constitutive_phenopowerlaw_Cslip_66 = 0.0_pReal + + allocate(constitutive_phenopowerlaw_gdot0_slip(maxNinstance)) ; constitutive_phenopowerlaw_gdot0_slip = 0.0_pReal + allocate(constitutive_phenopowerlaw_n_slip(maxNinstance)) ; constitutive_phenopowerlaw_n_slip = 0.0_pReal + allocate(constitutive_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,& + maxNinstance)) ; constitutive_phenopowerlaw_tau0_slip = 0.0_pReal + allocate(constitutive_phenopowerlaw_tausat_slip(lattice_maxNslipFamily,& + maxNinstance)) ; constitutive_phenopowerlaw_tausat_slip = 0.0_pReal + + allocate(constitutive_phenopowerlaw_gdot0_twin(maxNinstance)) ; constitutive_phenopowerlaw_gdot0_twin = 0.0_pReal + allocate(constitutive_phenopowerlaw_n_twin(maxNinstance)) ; constitutive_phenopowerlaw_n_twin = 0.0_pReal + allocate(constitutive_phenopowerlaw_tau0_twin(lattice_maxNtwinFamily,& + maxNinstance)) ; constitutive_phenopowerlaw_tau0_twin = 0.0_pReal + + allocate(constitutive_phenopowerlaw_spr(maxNinstance)) ; constitutive_phenopowerlaw_spr = 0.0_pReal + allocate(constitutive_phenopowerlaw_twinB(maxNinstance)) ; constitutive_phenopowerlaw_twinB = 0.0_pReal + allocate(constitutive_phenopowerlaw_twinC(maxNinstance)) ; constitutive_phenopowerlaw_twinC = 0.0_pReal + allocate(constitutive_phenopowerlaw_twinD(maxNinstance)) ; constitutive_phenopowerlaw_twinD = 0.0_pReal + allocate(constitutive_phenopowerlaw_twinE(maxNinstance)) ; constitutive_phenopowerlaw_twinE = 0.0_pReal + + allocate(constitutive_phenopowerlaw_h0_slipslip(maxNinstance)) ; constitutive_phenopowerlaw_h0_slipslip = 0.0_pReal + allocate(constitutive_phenopowerlaw_h0_sliptwin(maxNinstance)) ; constitutive_phenopowerlaw_h0_sliptwin = 0.0_pReal + allocate(constitutive_phenopowerlaw_h0_twinslip(maxNinstance)) ; constitutive_phenopowerlaw_h0_twinslip = 0.0_pReal + allocate(constitutive_phenopowerlaw_h0_twintwin(maxNinstance)) ; constitutive_phenopowerlaw_h0_twintwin = 0.0_pReal + + allocate(constitutive_phenopowerlaw_interaction_slipslip(lattice_maxNinteraction,& + maxNinstance)) ; constitutive_phenopowerlaw_interaction_slipslip = 0.0_pReal + allocate(constitutive_phenopowerlaw_interaction_sliptwin(lattice_maxNinteraction,& + maxNinstance)) ; constitutive_phenopowerlaw_interaction_sliptwin = 0.0_pReal + allocate(constitutive_phenopowerlaw_interaction_twinslip(lattice_maxNinteraction,& + maxNinstance)) ; constitutive_phenopowerlaw_interaction_twinslip = 0.0_pReal + allocate(constitutive_phenopowerlaw_interaction_twintwin(lattice_maxNinteraction,& + maxNinstance)) ; constitutive_phenopowerlaw_interaction_twintwin = 0.0_pReal + + rewind(file) + line = '' + section = 0 + do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to + read(file,'(a1024)',END=100) line + enddo + + do ! read thru sections of phase part + read(file,'(a1024)',END=100) line + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') exit ! stop at next part + if (IO_getTag(line,'[',']') /= '') then ! next section + section = section + 1 + output = 0 ! reset output counter + endif + if (section > 0 .and. phase_constitution(section) == constitutive_phenopowerlaw_label) then ! one of my sections + i = phase_constitutionInstance(section) ! which instance of my constitution is present phase + positions = IO_stringPos(line,maxNchunks) + tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key + select case(tag) + case ('(output)') + output = output + 1 + constitutive_phenopowerlaw_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) + case ('lattice_structure') + constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2)) + case ('covera_ratio') + constitutive_phenopowerlaw_CoverA(i) = IO_floatValue(line,positions,2) + case ('c11') + constitutive_phenopowerlaw_C11(i) = IO_floatValue(line,positions,2) + case ('c12') + constitutive_phenopowerlaw_C12(i) = IO_floatValue(line,positions,2) + case ('c13') + constitutive_phenopowerlaw_C13(i) = IO_floatValue(line,positions,2) + case ('c33') + constitutive_phenopowerlaw_C33(i) = IO_floatValue(line,positions,2) + case ('c44') + constitutive_phenopowerlaw_C44(i) = IO_floatValue(line,positions,2) + case ('nslip') + forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_Nslip(j,i) = IO_intValue(line,positions,1+j) + case ('gdot0_slip') + constitutive_phenopowerlaw_gdot0_slip(i) = IO_floatValue(line,positions,2) + case ('n_slip') + constitutive_phenopowerlaw_n_slip(i) = IO_floatValue(line,positions,2) + case ('tau0_slip') + forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_tau0_slip(j,i) = IO_floatValue(line,positions,1+j) + case ('tausat_slip') + forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_tausat_slip(j,i) = IO_floatValue(line,positions,1+j) + case ('ntwin') + forall (j = 1:lattice_maxNtwinFamily) constitutive_phenopowerlaw_Ntwin(j,i) = IO_intValue(line,positions,1+j) + case ('gdot0_twin') + constitutive_phenopowerlaw_gdot0_twin(i) = IO_floatValue(line,positions,2) + case ('n_twin') + constitutive_phenopowerlaw_n_twin(i) = IO_floatValue(line,positions,2) + case ('tau0_twin') + forall (j = 1:lattice_maxNtwinFamily) constitutive_phenopowerlaw_tau0_twin(j,i) = IO_floatValue(line,positions,1+j) + case ('s_pr') + constitutive_phenopowerlaw_spr(i) = IO_floatValue(line,positions,2) + case ('twin_b') + constitutive_phenopowerlaw_twinB(i) = IO_floatValue(line,positions,2) + case ('twin_c') + constitutive_phenopowerlaw_twinC(i) = IO_floatValue(line,positions,2) + case ('twin_d') + constitutive_phenopowerlaw_twinD(i) = IO_floatValue(line,positions,2) + case ('twin_e') + constitutive_phenopowerlaw_twinE(i) = IO_floatValue(line,positions,2) + case ('h0_slipslip') + constitutive_phenopowerlaw_h0_slipslip(i) = IO_floatValue(line,positions,2) + case ('h0_sliptwin') + constitutive_phenopowerlaw_h0_sliptwin(i) = IO_floatValue(line,positions,2) + case ('h0_twinslip') + constitutive_phenopowerlaw_h0_twinslip(i) = IO_floatValue(line,positions,2) + case ('h0_twintwin') + constitutive_phenopowerlaw_h0_twintwin(i) = IO_floatValue(line,positions,2) + case ('interaction_slipslip') + forall (j = 1:lattice_maxNinteraction) & + constitutive_phenopowerlaw_interaction_slipslip(j,i) = IO_floatValue(line,positions,1+j) + case ('interaction_sliptwin') + forall (j = 1:lattice_maxNinteraction) & + constitutive_phenopowerlaw_interaction_sliptwin(j,i) = IO_floatValue(line,positions,1+j) + case ('interaction_twinslip') + forall (j = 1:lattice_maxNinteraction) & + constitutive_phenopowerlaw_interaction_twinslip(j,i) = IO_floatValue(line,positions,1+j) + case ('interaction_twintwin') + forall (j = 1:lattice_maxNinteraction) & + constitutive_phenopowerlaw_interaction_twintwin(j,i) = IO_floatValue(line,positions,1+j) + end select + endif + enddo + +100 do i = 1,maxNinstance + constitutive_phenopowerlaw_structure(i) = lattice_initializeStructure(constitutive_phenopowerlaw_structureName(i), & ! get structure + constitutive_phenopowerlaw_CoverA(i)) + constitutive_phenopowerlaw_Nslip(:,i) = min(lattice_NslipSystem(:,constitutive_phenopowerlaw_structure(i)),& ! limit active slip systems per family to max + constitutive_phenopowerlaw_Nslip(:,i)) + constitutive_phenopowerlaw_Ntwin(:,i) = min(lattice_NtwinSystem(:,constitutive_phenopowerlaw_structure(i)),& ! limit active twin systems per family to max + constitutive_phenopowerlaw_Ntwin(:,i)) + constitutive_phenopowerlaw_totalNslip(i) = sum(constitutive_phenopowerlaw_Nslip(:,i)) ! how many slip systems altogether + constitutive_phenopowerlaw_totalNtwin(i) = sum(constitutive_phenopowerlaw_Ntwin(:,i)) ! how many twin systems altogether + + if (constitutive_phenopowerlaw_structure(i) < 1 .or. & ! sanity checks + constitutive_phenopowerlaw_structure(i) > 3) call IO_error(205) + + if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. & + constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(210) + if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211) + if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(212) + if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. & + constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(213) + + if (any(constitutive_phenopowerlaw_tau0_twin(:,i) < 0.0_pReal .and. & + constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(210) + if (constitutive_phenopowerlaw_gdot0_twin(i) <= 0.0_pReal) call IO_error(211) + if (constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal) call IO_error(212) + + enddo + + allocate(constitutive_phenopowerlaw_hardeningMatrix_slipslip(maxval(constitutive_phenopowerlaw_totalNslip),& + maxval(constitutive_phenopowerlaw_totalNslip),& + maxNinstance)) + allocate(constitutive_phenopowerlaw_hardeningMatrix_sliptwin(maxval(constitutive_phenopowerlaw_totalNslip),& + maxval(constitutive_phenopowerlaw_totalNtwin),& + maxNinstance)) + allocate(constitutive_phenopowerlaw_hardeningMatrix_twinslip(maxval(constitutive_phenopowerlaw_totalNtwin),& + maxval(constitutive_phenopowerlaw_totalNslip),& + maxNinstance)) + allocate(constitutive_phenopowerlaw_hardeningMatrix_twintwin(maxval(constitutive_phenopowerlaw_totalNtwin),& + maxval(constitutive_phenopowerlaw_totalNtwin),& + maxNinstance)) + + do i = 1,maxNinstance + do j = 1,maxval(phase_Noutput) + select case(constitutive_phenopowerlaw_output(j,i)) + case('resistance_slip') + mySize = constitutive_phenopowerlaw_totalNslip(i) + case('shearrate_slip') + mySize = constitutive_phenopowerlaw_totalNslip(i) + case('resolvedstress_slip') + mySize = constitutive_phenopowerlaw_totalNslip(i) + case('totalshear') + mySize = 1_pInt + case('resistance_twin') + mySize = constitutive_phenopowerlaw_totalNtwin(i) + case('shearrate_twin') + mySize = constitutive_phenopowerlaw_totalNtwin(i) + case('resolvedstress_twin') + mySize = constitutive_phenopowerlaw_totalNtwin(i) + case('totalvolfrac') + mySize = 1_pInt + case default + mySize = 0_pInt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + constitutive_phenopowerlaw_sizePostResult(j,i) = mySize + constitutive_phenopowerlaw_sizePostResults(i) = & + constitutive_phenopowerlaw_sizePostResults(i) + mySize + endif + enddo + + constitutive_phenopowerlaw_sizeDotState(i) = constitutive_phenopowerlaw_totalNslip(i)+ & + constitutive_phenopowerlaw_totalNtwin(i)+ 2 ! s_slip, s_twin, sum(gamma), sum(f) + 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) + 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 + 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) + constitutive_phenopowerlaw_Cslip_66(1,2,i) = constitutive_phenopowerlaw_C12(i) + constitutive_phenopowerlaw_Cslip_66(2,1,i) = constitutive_phenopowerlaw_C12(i) + constitutive_phenopowerlaw_Cslip_66(1,3,i) = constitutive_phenopowerlaw_C13(i) + constitutive_phenopowerlaw_Cslip_66(3,1,i) = constitutive_phenopowerlaw_C13(i) + constitutive_phenopowerlaw_Cslip_66(2,3,i) = constitutive_phenopowerlaw_C13(i) + constitutive_phenopowerlaw_Cslip_66(3,2,i) = constitutive_phenopowerlaw_C13(i) + constitutive_phenopowerlaw_Cslip_66(4,4,i) = constitutive_phenopowerlaw_C44(i) + constitutive_phenopowerlaw_Cslip_66(5,5,i) = constitutive_phenopowerlaw_C44(i) + constitutive_phenopowerlaw_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_phenopowerlaw_C11(i)- & + constitutive_phenopowerlaw_C12(i)) + end select + constitutive_phenopowerlaw_Cslip_66(:,:,i) = & + math_Mandel3333to66(math_Voigt66to3333(constitutive_phenopowerlaw_Cslip_66(:,:,i))) + + do j = 1,lattice_maxNslipFamily + do k = 1,constitutive_phenopowerlaw_Nslip(j,i) + do l = 1,lattice_maxNslipFamily + do m = 1,constitutive_phenopowerlaw_Nslip(l,i) + constitutive_phenopowerlaw_hardeningMatrix_slipslip(sum(constitutive_phenopowerlaw_Nslip(1:j-1,i))+k,& + sum(constitutive_phenopowerlaw_Nslip(1:l-1,i))+m, i) = & + constitutive_phenopowerlaw_interaction_slipslip(lattice_interactionSlipSlip( & + sum(lattice_NslipSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, & + sum(lattice_NslipSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, & + constitutive_phenopowerlaw_structure(i)), i ) + enddo; enddo; enddo; enddo + + do j = 1,lattice_maxNslipFamily + do k = 1,constitutive_phenopowerlaw_Nslip(j,i) + do l = 1,lattice_maxNtwinFamily + do m = 1,constitutive_phenopowerlaw_Ntwin(l,i) + constitutive_phenopowerlaw_hardeningMatrix_sliptwin(sum(constitutive_phenopowerlaw_Nslip(1:j-1,i))+k,& + sum(constitutive_phenopowerlaw_Ntwin(1:l-1,i))+m, i) = & + constitutive_phenopowerlaw_interaction_sliptwin(lattice_interactionSlipTwin( & + sum(lattice_NslipSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, & + sum(lattice_NtwinSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, & + constitutive_phenopowerlaw_structure(i)) ,i ) + enddo; enddo; enddo; enddo + + do j = 1,lattice_maxNtwinFamily + do k = 1,constitutive_phenopowerlaw_Ntwin(j,i) + do l = 1,lattice_maxNslipFamily + do m = 1,constitutive_phenopowerlaw_Nslip(l,i) + constitutive_phenopowerlaw_hardeningMatrix_twinslip(sum(constitutive_phenopowerlaw_Ntwin(1:j-1,i))+k,& + sum(constitutive_phenopowerlaw_Nslip(1:l-1,i))+m, i) = & + constitutive_phenopowerlaw_interaction_twinslip(lattice_interactionTwinSlip( & + sum(lattice_NtwinSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, & + sum(lattice_NslipSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, & + constitutive_phenopowerlaw_structure(i)), i ) + enddo; enddo; enddo; enddo + + do j = 1,lattice_maxNtwinFamily + do k = 1,constitutive_phenopowerlaw_Ntwin(j,i) + do l = 1,lattice_maxNtwinFamily + do m = 1,constitutive_phenopowerlaw_Ntwin(l,i) + constitutive_phenopowerlaw_hardeningMatrix_twintwin(sum(constitutive_phenopowerlaw_Ntwin(1:j-1,i))+k,& + sum(constitutive_phenopowerlaw_Ntwin(1:l-1,i))+m, i) = & + constitutive_phenopowerlaw_interaction_twintwin(lattice_interactionTwinTwin( & + sum(lattice_NtwinSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, & + sum(lattice_NtwinSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, & + constitutive_phenopowerlaw_structure(i)), i ) + enddo; enddo; enddo; enddo + + enddo + + return + +end subroutine + + +function constitutive_phenopowerlaw_stateInit(myInstance) +!********************************************************************* +!* initial microstructural state * +!********************************************************************* + use prec, only: pReal,pInt + use debug, only: debugger + use lattice, only: lattice_maxNslipFamily, lattice_maxNtwinFamily + implicit none + +!* Definition of variables + integer(pInt), intent(in) :: myInstance + integer(pInt) i + real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(myInstance) + & + constitutive_phenopowerlaw_totalNtwin(myInstance) + 2) :: constitutive_phenopowerlaw_stateInit + + constitutive_phenopowerlaw_stateInit = 0.0_pReal + + do i = 1,lattice_maxNslipFamily + constitutive_phenopowerlaw_stateInit(1+& + sum(constitutive_phenopowerlaw_Nslip(1:i-1,myInstance)) : & + sum(constitutive_phenopowerlaw_Nslip(1:i ,myInstance))) = & + constitutive_phenopowerlaw_tau0_slip(i,myInstance) + enddo + + do i = 1,lattice_maxNtwinFamily + constitutive_phenopowerlaw_stateInit(1+sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+& + sum(constitutive_phenopowerlaw_Ntwin(1:i-1,myInstance)) : & + sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+& + sum(constitutive_phenopowerlaw_Ntwin(1:i ,myInstance))) = & + constitutive_phenopowerlaw_tau0_twin(i,myInstance) + enddo + return + +end function + + +function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el) +!********************************************************************* +!* homogenized elacticity matrix * +!* INPUT: * +!* - state : state variables * +!* - ipc : component-ID of current integration point * +!* - ip : current integration point * +!* - el : current element * +!********************************************************************* + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance + implicit none + +!* Definition of variables + integer(pInt), intent(in) :: ipc,ip,el + integer(pInt) matID + real(pReal), dimension(6,6) :: constitutive_phenopowerlaw_homogenizedC + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state + + matID = phase_constitutionInstance(material_phase(ipc,ip,el)) + constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(:,:,matID) + + return + +end function + + +subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el) +!********************************************************************* +!* calculate derived quantities from state (not used here) * +!* INPUT: * +!* - Tp : temperature * +!* - ipc : component-ID of current integration point * +!* - ip : current integration point * +!* - el : current element * +!********************************************************************* + use prec, only: pReal,pInt,p_vec + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance + implicit none + +!* Definition of variables + integer(pInt) ipc,ip,el, matID,nSlip,nTwin + real(pReal) Temperature + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state + + matID = phase_constitutionInstance(material_phase(ipc,ip,el)) + +end subroutine + + +subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el) +!********************************************************************* +!* plastic velocity gradient and its tangent * +!* INPUT: * +!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * +!* - ipc : component-ID at current integration point * +!* - ip : current integration point * +!* - el : current element * +!* OUTPUT: * +!* - Lp : plastic velocity gradient * +!* - dLp_dTstar : derivative of Lp (4th-rank tensor) * +!********************************************************************* + use prec, only: pReal,pInt,p_vec + use debug, only: debugger + use math, only: math_Plain3333to99 + use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & + lattice_NslipSystem,lattice_NtwinSystem + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance + + implicit none + +!* Definition of variables + integer(pInt) ipc,ip,el + integer(pInt) matID,nSlip,nTwin,f,i,j,k,l,m,n, structID,index_Gamma,index_F,index_myFamily + real(pReal) Temperature + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(3,3) :: Lp + real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 + real(pReal), dimension(9,9) :: dLp_dTstar + real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip,dgdot_dtauslip,tau_slip + real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + gdot_twin,dgdot_dtautwin,tau_twin + + matID = phase_constitutionInstance(material_phase(ipc,ip,el)) + structID = constitutive_phenopowerlaw_structure(matID) + + nSlip = constitutive_phenopowerlaw_totalNslip(matID) + nTwin = constitutive_phenopowerlaw_totalNtwin(matID) + + index_Gamma = nSlip + nTwin + 1 + index_F = nSlip + nTwin + 2 + + Lp = 0.0_pReal + dLp_dTstar3333 = 0.0_pReal + dLp_dTstar = 0.0_pReal + j = 0_pInt + do f = 1,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + j = j+1_pInt + +!* Calculation of Lp + + tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID)) + gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**& + constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j)) + Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F + gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,structID) + +!* Calculation of the tangent of Lp + + dgdot_dtauslip(j) = gdot_slip(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip(j) + 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) + & + dgdot_dtauslip(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* & + lattice_Sslip(m,n,index_myFamily+i,structID) + enddo + enddo + + j = 0_pInt + do f = 1,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + j = j+1_pInt + +!* Calculation of Lp + + tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID)) + gdot_twin(j) = constitutive_phenopowerlaw_gdot0_twin(matID)*& + (abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**& + constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) + Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,structID) + +!* Calculation of the tangent of Lp + + dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(matID)/tau_twin(j) + 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) + & + dgdot_dtautwin(j)*lattice_Stwin(k,l,index_myFamily+i,structID)* & + lattice_Stwin(m,n,index_myFamily+i,structID) + enddo + enddo + + dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) + + return +end subroutine + + +function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el) +!********************************************************************* +!* rate of change of microstructure * +!* INPUT: * +!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * +!* - ipc : component-ID at current integration point * +!* - ip : current integration point * +!* - el : current element * +!* OUTPUT: * +!* - constitutive_dotState : evolution of state variable * +!********************************************************************* + use prec, only: pReal,pInt,p_vec + use debug, only: debugger + use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & + lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance + implicit none + +!* Definition of variables + integer(pInt) ipc,ip,el + integer(pInt) matID,nSlip,nTwin,f,i,j,k, structID,index_Gamma,index_F,index_myFamily + real(pReal) Temperature,c_slipslip,c_sliptwin,c_twinslip,c_twintwin, ssat + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state + real(pReal), dimension(6) :: Tstar_v + real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip,tau_slip,h_slipslip,h_sliptwin,N_slipslip,N_twinslip + real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + gdot_twin,tau_twin,h_twinslip,h_twintwin,N_sliptwin,N_twintwin + real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + constitutive_phenopowerlaw_dotState + + matID = phase_constitutionInstance(material_phase(ipc,ip,el)) + structID = constitutive_phenopowerlaw_structure(matID) + + nSlip = constitutive_phenopowerlaw_totalNslip(matID) + nTwin = constitutive_phenopowerlaw_totalNtwin(matID) + + index_Gamma = nSlip + nTwin + 1 + index_F = nSlip + nTwin + 2 + + constitutive_phenopowerlaw_dotState = 0.0_pReal + +!-- system-independent (nonlinear) prefactors to M_xx matrices + + c_slipslip = constitutive_phenopowerlaw_h0_slipslip(matID)*& + (1.0_pReal + & + constitutive_phenopowerlaw_twinC(matID)*state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinB(matID)) + c_sliptwin = 0.0_pReal + c_twinslip = constitutive_phenopowerlaw_h0_twinslip(matID)*& + state(ipc,ip,el)%p(index_Gamma)**constitutive_phenopowerlaw_twinE(matID) + c_twintwin = constitutive_phenopowerlaw_h0_twintwin(matID)*& + state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(matID) + +!-- add system-dependent part and calculate dot gammas + + j = 0_pInt + do f = 1,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family + ssat = constitutive_phenopowerlaw_tausat_slip(f,matID) + & + constitutive_phenopowerlaw_spr(matID)*dsqrt(state(ipc,ip,el)%p(index_F)) + do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + j = j+1_pInt + h_slipslip(j) = c_slipslip*(1.0_pReal-state(ipc,ip,el)%p(j)/ssat) ! system-dependent prefactor for slip--slip interaction + + h_sliptwin(j) = c_sliptwin ! no system-dependent part + +!* Calculation of dot gamma + + tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID)) + gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**& + constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j)) + enddo + enddo + + j = 0_pInt + do f = 1,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + j = j+1_pInt + h_twinslip(j) = c_twinslip ! no system-dependent parts + h_twintwin(j) = c_twintwin + +!* Calculation of dot vol frac + + tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID)) + gdot_twin(j) = constitutive_phenopowerlaw_gdot0_twin(matID)*& + (abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**& + constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) + enddo + enddo + +!-- calculate the overall hardening based on above + + j = 0_pInt + do f = 1,lattice_maxNslipFamily ! loop over all slip families + do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + j = j+1_pInt + forall (k=1:nSlip) N_slipslip(k) = constitutive_phenopowerlaw_hardeningMatrix_slipslip(j,k,matID) * & + abs(gdot_slip(k)) ! dot gamma_slip + forall (k=1:nTwin) N_sliptwin(k) = constitutive_phenopowerlaw_hardeningMatrix_sliptwin(j,k,matID) * & + gdot_twin(k) ! dot gamma_twin + constitutive_phenopowerlaw_dotState(j) = h_slipslip(j)*sum(N_slipslip) + & ! evolution of slip resistance j + h_sliptwin(j)*sum(N_sliptwin) + constitutive_phenopowerlaw_dotState(index_Gamma) = constitutive_phenopowerlaw_dotState(index_Gamma) + & + abs(gdot_slip(j)) + enddo + enddo + + j = 0_pInt + do f = 1,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + j = j+1_pInt + forall (k=1:nSlip) N_twinslip(k) = constitutive_phenopowerlaw_hardeningMatrix_twinslip(j,k,matID) * & + gdot_slip(k) ! dot gamma_slip + forall (k=1:nTwin) N_twintwin(k) = constitutive_phenopowerlaw_hardeningMatrix_twintwin(j,k,matID) * & + gdot_twin(k) ! dot gamma_twin + constitutive_phenopowerlaw_dotState(j+nSlip) = h_twinslip(j)*sum(N_twinslip) + & ! evolution of twin resistance j + h_twintwin(j)*sum(N_twintwin) + constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + & + gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID) + enddo + enddo + + return + +end function + + +!**************************************************************** +!* calculates the rate of change of temperature * +!**************************************************************** +pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el) + + !*** variables and functions from other modules ***! + use prec, only: pReal,pInt,p_vec + use lattice, only: lattice_Sslip_v + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance + implicit none + + !*** input variables ***! + real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: Temperature + integer(pInt), intent(in):: ipc, & ! grain number + ip, & ! integration point number + el ! element number + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure + + !*** output variables ***! + real(pReal) constitutive_phenopowerlaw_dotTemperature ! rate of change of temparature + + ! calculate dotTemperature + constitutive_phenopowerlaw_dotTemperature = 0.0_pReal + + return +endfunction + + + +pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) +!********************************************************************* +!* return array of constitutive results * +!* INPUT: * +!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) * +!* - dt : current time increment * +!* - ipc : component-ID at current integration point * +!* - ip : current integration point * +!* - el : current element * +!********************************************************************* + use prec, only: pReal,pInt,p_vec + use lattice, only: lattice_Sslip_v,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & + lattice_NslipSystem,lattice_NtwinSystem + use mesh, only: mesh_NcpElems,mesh_maxNips + use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput + implicit none + +!* Definition of variables + integer(pInt), intent(in) :: ipc,ip,el + real(pReal), intent(in) :: dt,Temperature + real(pReal), dimension(6), intent(in) :: Tstar_v + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state + integer(pInt) matID,o,f,i,c,nSlip,nTwin,j, structID,index_Gamma,index_F,index_myFamily + real(pReal) tau + real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: & + constitutive_phenopowerlaw_postResults + + matID = phase_constitutionInstance(material_phase(ipc,ip,el)) + structID = constitutive_phenopowerlaw_structure(matID) + + nSlip = constitutive_phenopowerlaw_totalNslip(matID) + nTwin = constitutive_phenopowerlaw_totalNtwin(matID) + + index_Gamma = nSlip + nTwin + 1 + index_F = nSlip + nTwin + 2 + + constitutive_phenopowerlaw_postResults = 0.0_pReal + c = 0_pInt + + do o = 1,phase_Noutput(material_phase(ipc,ip,el)) + select case(constitutive_phenopowerlaw_output(o,matID)) + case ('resistance_slip') + constitutive_phenopowerlaw_postResults(c+1:c+nSlip) = state(ipc,ip,el)%p(1:nSlip) + c = c + nSlip + + case ('shearrate_slip') + j = 0_pInt + do f = 1,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + j = j + 1_pInt + tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID)) + constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(matID)*& + (abs(tau)/state(ipc,ip,el)%p(j))**& + constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau) + enddo; enddo + c = c + nSlip + + case ('resolvedstress_slip') + j = 0_pInt + do f = 1,lattice_maxNslipFamily ! loop over all slip families + index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family + j = j + 1_pInt + constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID)) + enddo; enddo + c = c + nSlip + + case ('totalshear') + constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_Gamma) + c = c + 1 + + case ('resistance_twin') + constitutive_phenopowerlaw_postResults(c+1:c+nTwin) = state(ipc,ip,el)%p(1+nSlip:nTwin+nSlip) + c = c + nTwin + + case ('shearrate_twin') + j = 0_pInt + do f = 1,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + j = j + 1_pInt + tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID)) + constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_twin(matID)*& + (abs(tau)/state(ipc,ip,el)%p(j+nSlip))**& + constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau)) + enddo; enddo + c = c + nTwin + + case ('resolvedstress_twin') + j = 0_pInt + do f = 1,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family + do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family + j = j + 1_pInt + constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID)) + enddo; enddo + c = c + nTwin + + case ('totalvolfrac') + constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_F) + c = c + 1 + + end select + enddo + + return + +end function + +END MODULE diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 57474e901..2779a8d29 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -25,9 +25,9 @@ real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & crystallite_subdt, & ! substepped time increment of each grain crystallite_subFrac, & ! already calculated fraction of increment crystallite_subStep, & ! size of next integration step - crystallite_Temperature, & ! Temp of each grain - crystallite_partionedTemperature0, & ! Temp of each grain at start of homog inc - crystallite_subTemperature0 ! Temp of each grain at start of crystallite inc + crystallite_Temperature, & ! Temp of each grain + crystallite_partionedTemperature0, & ! Temp of each grain at start of homog inc + crystallite_subTemperature0 ! Temp of each grain at start of crystallite inc real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v, & ! current 2nd Piola-Kirchhoff stress vector (end of converged time step) crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partionedTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of homog inc @@ -55,7 +55,6 @@ logical, dimension (:,:,:), allocatable :: crystallite_localConstit crystallite_onTrack, & ! flag to indicate ongoing calculation crystallite_converged ! convergence flag - CONTAINS !******************************************************************** @@ -82,7 +81,7 @@ subroutine crystallite_init(Temperature) phase_localConstitution implicit none - !*** input variables ***! + !*** input variables ***! real(pReal) Temperature !*** output variables ***! @@ -95,40 +94,40 @@ subroutine crystallite_init(Temperature) iMax, & ! maximum number of integration points eMax, & ! maximum number of elements myNgrains - + gMax = homogenization_maxNgrains iMax = mesh_maxNips eMax = mesh_NcpElems - allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature - allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal + allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = Temperature + allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal - allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal + allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal - allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal - allocate(crystallite_partionedTemperature0(gMax,iMax,eMax)); crystallite_partionedTemperature0 = 0.0_pReal + allocate(crystallite_Tstar0_v(6,gMax,iMax,eMax)); crystallite_Tstar0_v = 0.0_pReal + allocate(crystallite_partionedTemperature0(gMax,iMax,eMax)); crystallite_partionedTemperature0 = 0.0_pReal allocate(crystallite_partionedF(3,3,gMax,iMax,eMax)); crystallite_partionedF = 0.0_pReal allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax)); crystallite_partionedFp0 = 0.0_pReal allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal - allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal - allocate(crystallite_subTemperature0(gMax,iMax,eMax)); crystallite_subTemperature0 = 0.0_pReal + allocate(crystallite_partionedTstar0_v(6,gMax,iMax,eMax)); crystallite_partionedTstar0_v = 0.0_pReal + allocate(crystallite_subTemperature0(gMax,iMax,eMax)); crystallite_subTemperature0 = 0.0_pReal allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal - allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal + allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal allocate(crystallite_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal - allocate(crystallite_localConstitution(gMax,iMax,eMax)); + allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true. allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false. allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. @@ -138,16 +137,15 @@ subroutine crystallite_init(Temperature) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element do g = 1,myNgrains - crystallite_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption + 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_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) crystallite_partionedF(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) crystallite_requested(g,i,e) = .true. - crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) - - enddo + crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) + enddo enddo enddo !$OMPEND PARALLEL DO @@ -161,20 +159,20 @@ subroutine crystallite_init(Temperature) write(6,*) '<<<+- crystallite init -+>>>' write(6,*) write(6,'(a32,x,7(i5,x))') 'crystallite_Nresults: ', crystallite_Nresults - write(6,'(a32,x,7(i5,x))') 'crystallite_Temperature ', shape(crystallite_Temperature) - write(6,'(a32,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe) + write(6,'(a32,x,7(i5,x))') 'crystallite_Temperature ', shape(crystallite_Temperature) + write(6,'(a32,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe) write(6,'(a32,x,7(i5,x))') 'crystallite_Fp: ', shape(crystallite_Fp) write(6,'(a32,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp) - write(6,'(a32,x,7(i5,x))') 'crystallite_F0: ', shape(crystallite_F0) + write(6,'(a32,x,7(i5,x))') 'crystallite_F0: ', shape(crystallite_F0) write(6,'(a32,x,7(i5,x))') 'crystallite_Fp0: ', shape(crystallite_Fp0) write(6,'(a32,x,7(i5,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) - write(6,'(a32,x,7(i5,x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0) + write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF) + write(6,'(a32,x,7(i5,x))') 'crystallite_partionedTemp0: ', shape(crystallite_partionedTemperature0) write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0) write(6,'(a32,x,7(i5,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0) write(6,'(a32,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0) - write(6,'(a32,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF) - write(6,'(a32,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) + write(6,'(a32,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF) + write(6,'(a32,x,7(i5,x))') 'crystallite_subTemperature0: ', shape(crystallite_subTemperature0) write(6,'(a32,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0) write(6,'(a32,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0) write(6,'(a32,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0) @@ -250,8 +248,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !*** output variables ***! !*** local variables ***! - real(pReal) myTemperature ! local copy of the temperature - real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient + real(pReal) myTemperature ! local copy of the temperature + real(pReal), dimension(3,3) :: invFp, & ! inverse of the plastic deformation gradient Fe_guess, & ! guess for elastic deformation gradient Tstar, & ! 2nd Piola-Kirchhoff stress tensor myF, & ! local copy of the deformation gradient @@ -274,16 +272,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco) converged ! flag indicating if iteration converged - ! ------ initialize to starting condition ------ - write (6,*) - write (6,*) 'Crystallite request from Materialpoint' - write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemp0 of 1 1 1' ,crystallite_partionedTemperature0(1,1,1) - write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1' ,crystallite_partionedF0(1:3,:,1,1,1) - write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1' ,crystallite_partionedFp0(1:3,:,1,1,1) - write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1' ,crystallite_partionedF(1:3,:,1,1,1) - write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1' ,crystallite_partionedLp0(1:3,:,1,1,1) +!$OMP CRITICAL (write2out) +! write (6,*) +! write (6,*) 'Crystallite request from Materialpoint' +! write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemp0 of 1 1 1' ,crystallite_partionedTemperature0(1,1,1) +! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1' ,crystallite_partionedF0(1:3,:,1,1,1) +! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1' ,crystallite_partionedFp0(1:3,:,1,1,1) +! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1' ,crystallite_partionedF(1:3,:,1,1,1) +! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1' ,crystallite_partionedLp0(1:3,:,1,1,1) +!$OMPEND CRITICAL (write2out) !$OMP PARALLEL DO @@ -292,7 +291,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains if (crystallite_requested(g,i,e)) then ! initialize restoration point of ... - crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature + crystallite_subTemperature0(g,i,e) = crystallite_partionedTemperature0(g,i,e) ! ...temperature constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad @@ -329,36 +328,37 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains + debugger = (e == 1 .and. i == 1 .and. g == 1) if (crystallite_converged(g,i,e)) then + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,'(a21,f10.8,a33,f10.8,a35)') 'winding forward from ', & + crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', & + crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent' + write(6,*) + !$OMPEND CRITICAL (write2out) + endif crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e)) if (crystallite_subStep(g,i,e) > subStepMin) then - crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... - crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ...def grad + crystallite_subTemperature0(g,i,e) = crystallite_Temperature(g,i,e) ! wind forward... + crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! ...def grad crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress endif - if (debugger) then - !$OMP CRITICAL (write2out) - write(6,'(a21,f6.4,a28,f6.4,a35)') 'winding forward from ', & - crystallite_subFrac(g,i,e)-crystallite_subStep(g,i,e),' to new crystallite_subfrac ', & - crystallite_subFrac(g,i,e),' in crystallite_stressAndItsTangent' - write(6,*) - !$OMPEND CRITICAL (write2out) - endif else crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore... - crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature - crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad + crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature + crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress if (debugger) then !$OMP CRITICAL (write2out) - write(6,'(a78,f6.4)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& - crystallite_subStep(g,i,e) + write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& + crystallite_subStep(g,i,e) write(6,*) !$OMPEND CRITICAL (write2out) endif @@ -389,8 +389,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! based on constitutive_subState0 ! results in constitutive_state - if (debugger) write(6,*) 'state integration started' - !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) @@ -399,8 +397,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) if ( crystallite_requested(g,i,e) & .and. crystallite_onTrack(g,i,e) & .and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites - crystallite_converged(g,i,e) = crystallite_updateState(g,i,e) ! use former evolution rate - crystallite_converged(g,i,e) = crystallite_updateTemperature(g,i,e) ! use former evolution rate + crystallite_converged(g,i,e) = crystallite_updateState(g,i,e) ! use former evolution rate + crystallite_converged(g,i,e) = crystallite_updateTemperature(g,i,e) ! use former evolution rate crystallite_converged(g,i,e) = .false. ! force at least one iteration step even if state already converged endif enddo @@ -431,6 +429,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains + debugger = (e == 1 .and. i == 1 .and. g == 1) if ( crystallite_requested(g,i,e) & .and. crystallite_onTrack(g,i,e) & .and. .not. crystallite_converged(g,i,e) ) & ! all undone crystallites @@ -451,11 +450,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains + debugger = (e == 1 .and. i == 1 .and. g == 1) if ( crystallite_requested(g,i,e) & .and. crystallite_onTrack(g,i,e) & .and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites - crystallite_converged(g,i,e) = (crystallite_updateState(g,i,e).AND.& - crystallite_updateTemperature(g,i,e)) + crystallite_converged(g,i,e) = (crystallite_updateState(g,i,e) .and. & + crystallite_updateTemperature(g,i,e)) + if (debugger) write (6,*) g,i,e,'converged after updState',crystallite_converged(g,i,e) if (crystallite_converged(g,i,e)) then !$OMP CRITICAL (distributionState) debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1 @@ -473,8 +474,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo ! crystallite convergence loop - if (debugger) write(6,*) 'state integration converged' - enddo ! cutback loop @@ -504,65 +503,69 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! --+>> stiffness calculation <<+-- if(updateJaco) then ! Jacobian required - if (debugger) write (6,*) 'Stiffness calculation started' - !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do g = 1,myNgrains +! debugger = (g == 1 .and. i == 1 .and. e == 1) if (crystallite_converged(g,i,e)) then ! grain converged in above iteration mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state... myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics - myFp = crystallite_Fp(:,:,g,i,e) - myTemperature = crystallite_Temperature(g,i,e) + myFp = crystallite_Fp(:,:,g,i,e) + myTemperature = crystallite_Temperature(g,i,e) myFe = crystallite_Fe(:,:,g,i,e) myLp = crystallite_Lp(:,:,g,i,e) myTstar_v = crystallite_Tstar_v(:,g,i,e) myP = crystallite_P(:,:,g,i,e) if (debugger) then + !$OMP CRITICAL (write2out) write (6,*) '#############' write (6,*) 'central solution' write (6,*) '#############' write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',myP(1:3,:)/1e6 write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:) write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:) - write (6,'(a,/,f12.4)') 'state of 1 1 1',myState/1e6 + write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6 + !$OMPEND CRITICAL (write2out) endif do k = 1,3 ! perturbation... do l = 1,3 ! ...components crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component if (debugger) then + !$OMP CRITICAL (write2out) write (6,*) '=============' write (6,'(i1,x,i1)') k,l write (6,*) '=============' write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) + !$OMPEND CRITICAL (write2out) endif onTrack = .true. converged = .false. NiterationState = 0_pInt do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) NiterationState = NiterationState + 1_pInt - if (debugger) write (6,'(a4,x,i6)') 'loop',NiterationState onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) - if (onTrack) converged = crystallite_updateState(g,i,e).AND.& ! update state + if (onTrack) converged = crystallite_updateState(g,i,e).AND.& ! update state crystallite_updateTemperature(g,i,e) ! update temperature if (debugger) then + !$OMP CRITICAL (write2out) write (6,*) '-------------' write (6,'(l,x,l)') onTrack,converged write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6 - write (6,'(a,/,f12.4)') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 - write (6,'(a,/,f12.4)') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 + write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 + write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 + !$OMPEND CRITICAL (write2out) endif enddo if (converged) & ! converged state warrants stiffness update crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_P(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state... - crystallite_Temperature(g,i,e)= myTemperature ! ... temperature - crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics + crystallite_Temperature(g,i,e)= myTemperature ! ... temperature + crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics crystallite_Fe(:,:,g,i,e) = myFe crystallite_Lp(:,:,g,i,e) = myLp crystallite_Tstar_v(:,g,i,e) = myTstar_v @@ -573,7 +576,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMPEND CRITICAL (out) enddo enddo - if (debugger) write (6,'(a,/,9(9(f12.4,x)/))') 'dPdF/GPa',crystallite_dPdF(:,:,:,:,1,1,1)/1e9 else ! grain did not converged crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback @@ -583,154 +585,151 @@ subroutine crystallite_stressAndItsTangent(updateJaco) enddo !$OMPEND PARALLEL DO - if (debugger) write (6,*) 'Stiffness calculation finished' - endif endsubroutine -!******************************************************************** -! update the internal state of the constitutive law -! and tell whether state has converged -!******************************************************************** - function crystallite_updateState(& - g,& ! grain number - i,& ! integration point number - e & ! element number - ) - - !*** variables and functions from other modules ***! - use prec, only: pReal, & - pInt, & - pLongInt - use numerics, only: rTol_crystalliteState - use constitutive, only: constitutive_dotState, & - constitutive_sizeDotState, & - constitutive_subState0, & - constitutive_state - use debug, only: debugger, & - debug_cumDotStateCalls, & - debug_cumDotStateTicks - - !*** input variables ***! - integer(pInt), intent(in):: e, & ! element index - i, & ! integration point index - g ! grain index - - !*** output variables ***! - logical crystallite_updateState ! flag indicating if integration suceeded - - !*** local variables ***! - real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure - integer(pInt) mySize - integer(pLongInt) tick, & - tock, & - tickrate, & - maxticks - - !*** global variables ***! - ! crystallite_Tstar_v - ! crystallite_subdt - ! crystallite_Temperature - - mySize = constitutive_sizeDotState(g,i,e) - - ! calculate the residuum - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - & - crystallite_subdt(g,i,e) * constitutive_dotState(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) - 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 - - ! if NaN occured then return without changing the state - if (any(residuum/=residuum)) then - crystallite_updateState = .false. ! indicate state update failed - if (debugger) write(6,*) '::: updateState encountered NaN' - return - endif - - ! update the microstructure - constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum - - ! setting flag to true if state is below relative Tolerance, otherwise set it to false - crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)), & - constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState - - if (debugger) write(6,'(a,/,f12.4)') 'updated state: ', constitutive_state(g,i,e)%p(1) - - return - - endfunction - - -!******************************************************************** -! update the temperature of the grain -! and tell whether it has converged -!******************************************************************** - function crystallite_updateTemperature(& - g,& ! grain number - i,& ! integration point number - e & ! element number - ) - - !*** variables and functions from other modules ***! - use prec, only: pReal, & - pInt, & - pLongInt - use numerics, only: rTol_crystalliteTemperature - use constitutive, only: constitutive_dotTemperature - use debug, only: debugger, & - debug_cumDotTemperatureCalls, & - debug_cumDotTemperatureTicks - - !*** input variables ***! - integer(pInt), intent(in):: e, & ! element index - i, & ! integration point index - g ! grain index - - !*** output variables ***! - logical crystallite_updateTemperature ! flag indicating if integration suceeded - - !*** local variables ***! - real(pReal) residuum ! residuum from evolution of temperature - integer(pLongInt) tick, & - tock, & - tickrate, & - maxticks - - ! calculate the residuum - call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) - & - crystallite_subdt(g,i,e) * constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) - call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) - debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt - debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick - if (tock < tick) debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + maxticks - - ! if NaN occured then return without changing the state - if (residuum/=residuum) then - crystallite_updateTemperature = .false. ! indicate update failed - if (debugger) write(6,*) '::: updateTemperature encountered NaN' - return - endif - - ! update the microstructure - crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum - - ! setting flag to true if state is below relative Tolerance, otherwise set it to false - crystallite_updateTemperature = abs(residuum/crystallite_Temperature(g,i,e)) < rTol_crystalliteTemperature - - if (debugger) write(6,'(a,/,f12.4)') 'updated temperature: ', crystallite_Temperature(g,i,e) - - return - - endfunction - - +!******************************************************************** +! update the internal state of the constitutive law +! and tell whether state has converged +!******************************************************************** + function crystallite_updateState(& + g,& ! grain number + i,& ! integration point number + e & ! element number + ) + + !*** variables and functions from other modules ***! + use prec, only: pReal, & + pInt, & + pLongInt + use numerics, only: rTol_crystalliteState + use constitutive, only: constitutive_dotState, & + constitutive_sizeDotState, & + constitutive_subState0, & + constitutive_state + use debug, only: debugger, & + debug_cumDotStateCalls, & + debug_cumDotStateTicks + + !*** input variables ***! + integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index + + !*** output variables ***! + logical crystallite_updateState ! flag indicating if integration suceeded + + !*** local variables ***! + real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum ! residuum from evolution of microstructure + integer(pInt) mySize + integer(pLongInt) tick, & + tock, & + tickrate, & + maxticks + + mySize = constitutive_sizeDotState(g,i,e) + + ! calculate the residuum + call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) + residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - & + crystallite_subdt(g,i,e) * constitutive_dotState(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) + 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 + + ! if NaN occured then return without changing the state + if (any(residuum/=residuum)) then + crystallite_updateState = .false. ! indicate state update failed + !$OMP CRITICAL (write2out) + write(6,*) '::: updateState encountered NaN',e,i,g + !$OMPEND CRITICAL (write2out) + return + endif + + ! update the microstructure + constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum + + ! setting flag to true if state is below relative tolerance, otherwise set it to false + crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)), & + constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState + + return + + endfunction + + +!******************************************************************** +! update the temperature of the grain +! and tell whether it has converged +!******************************************************************** + function crystallite_updateTemperature(& + g,& ! grain number + i,& ! integration point number + e & ! element number + ) + + !*** variables and functions from other modules ***! + use prec, only: pReal, & + pInt, & + pLongInt + use numerics, only: rTol_crystalliteTemperature + use constitutive, only: constitutive_dotTemperature + use debug, only: debugger, & + debug_cumDotTemperatureCalls, & + debug_cumDotTemperatureTicks + + !*** input variables ***! + integer(pInt), intent(in):: e, & ! element index + i, & ! integration point index + g ! grain index + + !*** output variables ***! + logical crystallite_updateTemperature ! flag indicating if integration suceeded + + !*** local variables ***! + real(pReal) residuum ! residuum from evolution of temperature + integer(pLongInt) tick, & + tock, & + tickrate, & + maxticks + + ! calculate the residuum + call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) + residuum = crystallite_Temperature(g,i,e) - crystallite_subTemperature0(g,i,e) - & + crystallite_subdt(g,i,e) * constitutive_dotTemperature(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) + call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) + debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt + debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick + if (tock < tick) debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + maxticks + + ! if NaN occured then return without changing the state + if (residuum/=residuum) then + crystallite_updateTemperature = .false. ! indicate update failed + !$OMP CRITICAL (write2out) + write(6,*) '::: updateTemperature encountered NaN',e,i,g + !$OMPEND CRITICAL (write2out) + return + endif + + ! update the microstructure + crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum + + ! setting flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false + if (crystallite_Temperature(g,i,e) /= 0.0_pReal) then + crystallite_updateTemperature = abs(residuum/crystallite_Temperature(g,i,e)) < rTol_crystalliteTemperature + else + crystallite_updateTemperature = .true. + endif + + return + + endfunction + + !*********************************************************************** !*** calculation of stress (P) with time integration *** @@ -824,16 +823,6 @@ endsubroutine tickrate, & maxticks - !*** global variables ***! - ! crystallite_subF - ! crystallite_subFp0 - ! crystallite_Tstar_v - ! crystallite_Lp - ! crystallite_subdt - ! crystallite_Temperature - - if (debugger) write(6,*) '::: integrateStress started' - ! be pessimistic crystallite_integrateStress = .false. @@ -847,7 +836,9 @@ endsubroutine ! inversion of Fp_current... invFp_current = math_inv3x3(Fp_current) if (all(invFp_current == 0.0_pReal)) then ! ... failed? - if (debugger) write(6,*) '::: integrateStress failed on invFp_current inversion' + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress failed on invFp_current inversion',e,i,g + !$OMPEND CRITICAL (write2out) return endif @@ -858,6 +849,7 @@ endsubroutine ! get elasticity tensor C_66 = constitutive_homogenizedC(g,i,e) +! if (debugger) write(6,'(a,/,6(6(f10.4,x)/))') 'elasticity',C_66(1:6,:)/1e9 C = math_Mandel66to3333(C_66) ! start LpLoop with no acceleration @@ -872,10 +864,7 @@ LpLoop: do NiterationStress = NiterationStress + 1 ! too many loops required ? - if (NiterationStress > nStress) then - if (debugger) write(6,*) '::: integrateStress exceeded nStress loopcount' - return - endif + if (NiterationStress > nStress) return B = math_I3 - crystallite_subdt(g,i,e)*Lpguess BT = transpose(B) @@ -910,7 +899,9 @@ LpLoop: do ! NaN occured at regular speed? if (any(residuum/=residuum) .and. leapfrog == 1.0) then - if (debugger) write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',e,i,g + !$OMPEND CRITICAL (write2out) return ! something went wrong at accelerated speed? @@ -941,7 +932,11 @@ LpLoop: do invdRdLp = 0.0_pReal call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR if (error) then - if (debugger) write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress + !$OMPEND CRITICAL (write2out) + endif return endif endif @@ -965,7 +960,11 @@ LpLoop: do invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det call math_invert3x3(invFp_new,Fp_new,det,error) if (error) then - if (debugger) write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress + if (debugger) then + !$OMP CRITICAL (write2out) + write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress + !$OMPEND CRITICAL (write2out) + endif return endif Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe @@ -985,10 +984,12 @@ LpLoop: do ! set return flag to true crystallite_integrateStress = .true. if (debugger) then + !$OMP CRITICAL (write2out) write(6,*) '::: integrateStress converged at iteration', NiterationStress write(6,*) - write(6,'(a,/,3(3(e15.7,x)/))') 'P of 1 1 1',crystallite_P(:,:,1,1,1) - write(6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 1 1',crystallite_Lp(:,:,1,1,1) + write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6 + write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e) + !$OMP CRITICAL (write2out) endif !$OMP CRITICAL (distributionStress) @@ -1043,11 +1044,7 @@ function crystallite_postResults(& real(pReal), dimension(3,3) :: U, & R logical error - - !*** global variables ***! - ! crystallite_Nresults - ! crystallite_Fe - + if (crystallite_Nresults >= 2) then crystallite_postResults(1) = material_phase(g,i,e) crystallite_postResults(2) = material_volume(g,i,e) diff --git a/code/debug.f90 b/code/debug.f90 index 95e71e507..462b86bdc 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -81,7 +81,7 @@ endsubroutine if (debug_cumLpCalls > 0_pInt) then call system_clock(count_rate=tickrate) 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,f12.3)') 'total CPU time/s :',dble(debug_cumLpTicks)/tickrate endif write(6,*) write(6,'(a33,x,i12)') 'total calls to dotState :',debug_cumDotStateCalls @@ -89,7 +89,7 @@ endsubroutine 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 + write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotStateTicks)/tickrate endif write(6,*) write(6,'(a33,x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls @@ -97,7 +97,7 @@ endsubroutine call system_clock(count_rate=tickrate) write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',& dble(debug_cumDotTemperatureTicks)/tickrate/1.0e-6_pReal/debug_cumDotTemperatureCalls - write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotTemperatureTicks + write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotTemperatureTicks)/tickrate endif integral = 0_pInt @@ -109,7 +109,7 @@ endsubroutine write(6,'(i25,i10)') i,debug_StressLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',sum(debug_StressLoopDistribution),integral + write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StressLoopDistribution) integral = 0_pInt write(6,*) @@ -120,7 +120,7 @@ endsubroutine write(6,'(i25,i10)') i,debug_StateLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',sum(debug_StateLoopDistribution),integral + write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StateLoopDistribution) integral = 0_pInt write(6,*) @@ -131,7 +131,7 @@ endsubroutine write(6,'(i25,i10)') i,debug_StiffnessStateLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',sum(debug_StiffnessStateLoopDistribution),integral + write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution) integral = 0_pInt write(6,*) @@ -142,7 +142,7 @@ endsubroutine write(6,'(i25,i10)') i,debug_CrystalliteLoopDistribution(i) endif enddo - write(6,'(a15,i10,i10)') ' total',sum(debug_CrystalliteLoopDistribution),integral + write(6,'(a15,i10,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution) write(6,*) endsubroutine diff --git a/code/homogenization.f90 b/code/homogenization.f90 index af939344b..70c316765 100644 --- a/code/homogenization.f90 +++ b/code/homogenization.f90 @@ -62,7 +62,7 @@ subroutine homogenization_init(Temperature) use homogenization_isostrain ! use homogenization_RGC - real(pReal) Temperature + real(pReal) Temperature integer(pInt), parameter :: fileunit = 200 integer(pInt) e,i,g,myInstance @@ -119,7 +119,7 @@ subroutine homogenization_init(Temperature) homogenization_maxSizeState = maxval(homogenization_sizeState) homogenization_maxSizePostResults = maxval(homogenization_sizePostResults) - allocate(materialpoint_results( 1+homogenization_maxSizePostResults + & + allocate(materialpoint_results( 1+ 1+homogenization_maxSizePostResults + & ! grain count, homogSize, homogResult homogenization_maxNgrains*(1+crystallite_Nresults+constitutive_maxSizePostResults), mesh_maxNips,mesh_NcpElems)) @@ -181,14 +181,14 @@ subroutine materialpoint_stressAndItsTangent(& use constitutive, only: constitutive_state0, & constitutive_partionedState0, & constitutive_state - use crystallite, only: crystallite_Temperature, & + use crystallite, only: crystallite_Temperature, & crystallite_F0, & crystallite_Fp0, & crystallite_Fp, & crystallite_Lp0, & crystallite_Lp, & crystallite_Tstar0_v, & - crystallite_Tstar_v, & + crystallite_Tstar_v, & crystallite_partionedTemperature0, & crystallite_partionedF0, & crystallite_partionedF, & @@ -210,8 +210,8 @@ subroutine materialpoint_stressAndItsTangent(& write (6,*) write (6,*) 'Material Point start' - write (6,'(a,/,(f12.7,x))') 'Temp0 of 8 1' ,materialpoint_Temperature(8,1) - write (6,'(a,/,3(3(f12.7,x)/))') 'F0 of 8 1',materialpoint_F0(1:3,:,8,1) + write (6,'(a,/,(f12.7,x))') 'Temp0 of 8 1' ,materialpoint_Temperature(8,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'F0 of 8 1',materialpoint_F0(1:3,:,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'F of 8 1',materialpoint_F(1:3,:,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Fp0 of 1 8 1',crystallite_Fp0(1:3,:,1,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 8 1',crystallite_Lp0(1:3,:,1,8,1) @@ -223,8 +223,8 @@ subroutine materialpoint_stressAndItsTangent(& ! initialize restoration points of grain... forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state0(g,i,e)%p ! ...microstructures - crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures - crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads + crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures + crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress @@ -263,8 +263,8 @@ subroutine materialpoint_stressAndItsTangent(& if (materialpoint_subStep(i,e) > subStepMin) then ! wind forward grain starting point of... - crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature(1:myNgrains,i,e) ! ...temperatures - crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads + crystallite_partionedTemperature0(1:myNgrains,i,e) = crystallite_Temperature(1:myNgrains,i,e) ! ...temperatures + crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedTstar0_v(:,1:myNgrains,i,e) = crystallite_Tstar_v(:,1:myNgrains,i,e) ! ...2nd PK stress @@ -281,8 +281,8 @@ subroutine materialpoint_stressAndItsTangent(& materialpoint_subStep(i,e) = 0.5_pReal * materialpoint_subStep(i,e) ! restore... - crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures - crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads + crystallite_Temperature(1:myNgrains,i,e) = crystallite_partionedTemperature0(1:myNgrains,i,e) ! ...temperatures + crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_Lp(:,:,1:myNgrains,i,e) = crystallite_partionedLp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_Tstar_v(:,1:myNgrains,i,e) = crystallite_partionedTstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures @@ -348,7 +348,7 @@ subroutine materialpoint_stressAndItsTangent(& do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed if ( materialpoint_requested(i,e) .and. & .not. materialpoint_doneAndHappy(1,i,e)) then - materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e) + materialpoint_doneAndHappy(:,i,e) = homogenization_updateState(i,e) materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(:,i,e)) ! converged if done and happy endif enddo @@ -365,14 +365,14 @@ subroutine materialpoint_stressAndItsTangent(& !$OMP PARALLEL DO do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed - call homogenization_averageStressAndItsTangent(i,e) + call homogenization_averageStressAndItsTangent(i,e) call homogenization_averageTemperature(i,e) enddo enddo !$OMP END PARALLEL DO write (6,*) 'Material Point finished' - write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1) + write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 1 1',crystallite_Lp(1:3,:,1,1,1) ! how to deal with stiffness? return @@ -401,9 +401,10 @@ subroutine materialpoint_postResults(dt) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed c = 0_pInt + materialpoint_results(c+1,i,e) = myNgrains; c = c+1_pInt ! tell number of grains at materialpoint d = homogenization_sizePostResults(i,e) materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of homogenization results - if (d > 0_pInt) then ! any homogenization results to mention? + if (d > 0_pInt) then ! any homogenization results to mention? materialpoint_results(c+1:c+d,i,e) = & ! tell homogenization results homogenization_postResults(i,e); c = c+d endif @@ -484,65 +485,65 @@ function homogenization_updateState(& endfunction -!******************************************************************** -! derive average stress and stiffness from constituent quantities -!******************************************************************** -subroutine homogenization_averageStressAndItsTangent(& - ip, & ! integration point - el & ! element - ) - use prec, only: pReal,pInt - use mesh, only: mesh_element - use material, only: homogenization_type, homogenization_maxNgrains - use crystallite, only: crystallite_P,crystallite_dPdF - - use homogenization_isostrain - implicit none - - integer(pInt), intent(in) :: ip,el - - select case(homogenization_type(mesh_element(3,el))) - case (homogenization_isostrain_label) - call homogenization_isostrain_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), & - materialpoint_dPdF(:,:,:,:,ip,el),& - crystallite_P(:,:,:,ip,el), & - crystallite_dPdF(:,:,:,:,:,ip,el), & - ip, & - el) - end select - - return - -endsubroutine - - -!******************************************************************** -! derive average stress and stiffness from constituent quantities -!******************************************************************** -subroutine homogenization_averageTemperature(& - ip, & ! integration point - el & ! element - ) - use prec, only: pReal,pInt - use mesh, only: mesh_element - use material, only: homogenization_type, homogenization_maxNgrains - use crystallite, only: crystallite_Temperature - - use homogenization_isostrain - implicit none - - integer(pInt), intent(in) :: ip,el - - select case(homogenization_type(mesh_element(3,el))) - case (homogenization_isostrain_label) - materialpoint_Temperature(ip,el) = homogenization_isostrain_averageTemperature(crystallite_Temperature(:,ip,el), ip, el) - end select - - return - -endsubroutine - - +!******************************************************************** +! derive average stress and stiffness from constituent quantities +!******************************************************************** +subroutine homogenization_averageStressAndItsTangent(& + ip, & ! integration point + el & ! element + ) + use prec, only: pReal,pInt + use mesh, only: mesh_element + use material, only: homogenization_type, homogenization_maxNgrains + use crystallite, only: crystallite_P,crystallite_dPdF + + use homogenization_isostrain + implicit none + + integer(pInt), intent(in) :: ip,el + + select case(homogenization_type(mesh_element(3,el))) + case (homogenization_isostrain_label) + call homogenization_isostrain_averageStressAndItsTangent( materialpoint_P(:,:,ip,el), & + materialpoint_dPdF(:,:,:,:,ip,el),& + crystallite_P(:,:,:,ip,el), & + crystallite_dPdF(:,:,:,:,:,ip,el), & + ip, & + el) + end select + + return + +endsubroutine + + +!******************************************************************** +! derive average stress and stiffness from constituent quantities +!******************************************************************** +subroutine homogenization_averageTemperature(& + ip, & ! integration point + el & ! element + ) + use prec, only: pReal,pInt + use mesh, only: mesh_element + use material, only: homogenization_type, homogenization_maxNgrains + use crystallite, only: crystallite_Temperature + + use homogenization_isostrain + implicit none + + integer(pInt), intent(in) :: ip,el + + select case(homogenization_type(mesh_element(3,el))) + case (homogenization_isostrain_label) + materialpoint_Temperature(ip,el) = homogenization_isostrain_averageTemperature(crystallite_Temperature(:,ip,el), ip, el) + end select + + return + +endsubroutine + + !******************************************************************** ! return array of homogenization results for post file inclusion ! call only, if homogenization_sizePostResults(ip,el) > 0 !! diff --git a/code/lattice.f90 b/code/lattice.f90 index 26ff927b3..56640ab92 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -20,41 +20,52 @@ implicit none integer(pInt) lattice_Nhexagonal, & ! # of hexagonal lattice structure (from tag CoverA_ratio) lattice_Nstructure ! # of lattice structures (1: fcc,2: bcc,3+: hexagonal) +integer(pInt), parameter :: lattice_maxNslipFamily = 4 ! max # of slip system families over lattice structures +integer(pInt), parameter :: lattice_maxNtwinFamily = 4 ! max # of twin system families over lattice structures integer(pInt), parameter :: lattice_maxNslip = 48 ! max # of slip systems over lattice structures 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), pointer, dimension(:,:) :: interactionSlipSlip, & interactionSlipTwin, & - interactionTwinTwin, & - interactionTwinSlip + interactionTwinSlip, & + interactionTwinTwin -! Schmid matrices, normal, shear direction and nxd of slip systems +! Schmid matrices, normal, shear direction and d x n of slip systems real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Sslip real(pReal), allocatable, dimension(:,:,:) :: lattice_Sslip_v -real(pReal), allocatable, dimension(:,:,:) :: lattice_sn -real(pReal), allocatable, dimension(:,:,:) :: lattice_sd -real(pReal), allocatable, dimension(:,:,:) :: lattice_st +real(pReal), allocatable, dimension(:,:,:) :: lattice_sn, & + lattice_sd, & + lattice_st -! Rotation and Schmid matrices, normal, shear direction and nxd of twin systems +! rotation and Schmid matrices, normal, shear direction and d x n of twin systems real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Qtwin real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Stwin real(pReal), allocatable, dimension(:,:,:) :: lattice_Stwin_v -real(pReal), allocatable, dimension(:,:,:) :: lattice_tn -real(pReal), allocatable, dimension(:,:,:) :: lattice_td -real(pReal), allocatable, dimension(:,:,:) :: lattice_tt +real(pReal), allocatable, dimension(:,:,:) :: lattice_tn, & + lattice_td, & + lattice_tt + +! characteristic twin shear real(pReal), allocatable, dimension(:,:) :: lattice_shearTwin -integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip -integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipTwin -integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinTwin -integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip +! number of slip and twin systems in each family +integer(pInt), allocatable, dimension(:,:) :: lattice_NslipSystem, & + lattice_NtwinSystem +! interaction type of slip and twin systems among each other +integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, & + lattice_interactionSlipTwin, & + lattice_interactionTwinSlip, & + lattice_interactionTwinTwin !============================== fcc (1) ================================= - integer(pInt), parameter :: lattice_fcc_Nslip = 12_pInt - integer(pInt), parameter :: lattice_fcc_Ntwin = 12_pInt + integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_fcc_NslipSystem = (/12, 0, 0, 0/) + integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_fcc_NtwinSystem = (/12, 0, 0, 0/) + integer(pInt), parameter :: lattice_fcc_Nslip = 12 ! sum(lattice_fcc_NslipSystem) + integer(pInt), parameter :: lattice_fcc_Ntwin = 12 ! sum(lattice_fcc_NtwinSystem) integer(pInt) :: lattice_fcc_Nstructure = 0_pInt real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter :: lattice_fcc_systemSlip = & @@ -161,8 +172,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip !============================== bcc (2) ================================= - integer(pInt), parameter :: lattice_bcc_Nslip = 48_pInt - integer(pInt), parameter :: lattice_bcc_Ntwin = 12_pInt + integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_bcc_NslipSystem = (/12,12,24, 0/) + integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_bcc_NtwinSystem = (/12, 0, 0, 0/) + integer(pInt), parameter :: lattice_bcc_Nslip = 48 ! sum(lattice_bcc_NslipSystem) + integer(pInt), parameter :: lattice_bcc_Ntwin = 12 ! sum(lattice_bcc_NtwinSystem) integer(pInt) :: lattice_bcc_Nstructure = 0_pInt real(pReal), dimension(3+3,lattice_bcc_Nslip), parameter :: lattice_bcc_systemSlip = & @@ -352,8 +365,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip !============================== hex (3+) ================================= - integer(pInt), parameter :: lattice_hex_Nslip = 24_pInt - integer(pInt), parameter :: lattice_hex_Ntwin = 24_pInt + integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_hex_NslipSystem = (/ 3, 3, 6,12/) + integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_hex_NtwinSystem = (/ 6, 6, 6, 6/) + integer(pInt), parameter :: lattice_hex_Nslip = 24 ! sum(lattice_hex_NslipSystem) + integer(pInt), parameter :: lattice_hex_Ntwin = 24 ! sum(lattice_hex_NtwinSystem) integer(pInt) :: lattice_hex_Nstructure = 0_pInt !* sorted by YJ.Ro and Philip @@ -391,29 +406,29 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter :: lattice_hex_systemTwin = & reshape((/& - 0, 1, -1, 1, 0, -1, 1, 2, & ! <1011>{1012} Twin: shear 0.169 -1.26 tensile + 0, 1, -1, 1, 0, -1, 1, 2, & ! <1011>{1012} Twin: shear 0.169 -1.26 compression -1, 1, 0, 1, 1, -1, 0, 2, & -1, 0, 1, 1, 1, 0, -1, 2, & 0, -1, 1, 1, 0, 1, -1, 2, & 1, -1, 0, 1, -1, 1, 0, 2, & 1, 0, -1, 1, -1, 0, 1, 2, & - 2, -1, -1, -3, 2, -1, -1, 2, & ! <211-2>{2112} Twin: shear 0.224 1.19 compressive + 2, -1, -1, -3, 2, -1, -1, 2, & ! <211-2>{2112} Twin: shear 0.224 1.19 tension 1, 1, -2, -3, 1, 1, -2, 2, & -1, 2, -1, -3, -1, 2, -1, 2, & -2, 1, 1, -3, -2, 1, 1, 2, & -1, -1, 2, -3, -1, -1, 2, 2, & 1, -2, 1, -3, 1, -2, 1, 2, & - -2, 1, 1, 6, 2, -1, -1, 1, & ! <211-6>{2111} Twin: shear 0.628 -0.39 tensile + -2, 1, 1, 6, 2, -1, -1, 1, & ! <211-6>{2111} Twin: shear 0.628 -0.39 compression -1, -1, 2, 6, 1, 1, -2, 1, & 1, -2, 1, 6, -1, 2, -1, 1, & 2, -1, -1, 6, -2, 1, 1, 1, & 1, 1, -2, 6, -1, -1, 2, 1, & -1, 2, -1, 6, 1, -2, 1, 1, & - 0, -1, 1, -2, 0, -1, 1, 1, & ! <101-2>{1011} Twin: shear 0.103 1.09 compressive - 1, -1, 0, -2, 1, -1, 0, 1, & - 1, 0, -1, -2, 1, 0, -1, 1, & - 0, 1, -1, -2, 0, 1, -1, 1, & + 1, 0, -1, -2, 1, 0, -1, 1, & ! <101-2>{1011} Twin: shear 0.103 1.09 tension -1, 0, 1, -2, -1, 0, 1, 1, & + 0, 1, -1, -2, 0, 1, -1, 1, & + 0, -1, 1, -2, 0, -1, 1, 1, & + 1, -1, 0, -2, 1, -1, 0, 1, & -1, 1, 0, -2, -1, 1, 0, 1 & /),(/4+4,lattice_hex_Ntwin/)) !* Sort? Numbering of twin system follows Prof. Tom Bieler's scheme (to be consistent with his work); but numbering in data was restarted from 1 & @@ -446,131 +461,124 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip /),(/lattice_hex_Ntwin/)) !* four different interaction type matrix - !* 1. slip-slip interaction - !- wihtin the same slip familiy, self or latent hardening; 1 (self), 2 ~ 5 (latent for each slip; bas, pri, pyr_a, pyr_c+a in sequence) - !- interaction between different slip families; 6 ~ 11 - !* 2. slip-twin interaction: current scheme -- indirect effect of twin volume - !- all 1 - !* 3. twin-twin interaction - !- wihtin the same slip familiy, self or latent hardening; 1 (self), 2 ~ 5 - !- interaction between different slip families; 6 ~ 11 - !* 4. twin-slip interaction - !- T1 interacting slip; 1 ~ 4, C1 interacting slip; 5 ~ 8 - !- T2 interacting slip; 9 ~ 12, C1 interacting slip; 13 ~ 16 + !* 1. slip-slip interaction - 20 types + !* 2. slip-twin interaction - 16 types + !* 3. twin-twin interaction - 20 types + !* 4. twin-slip interaction - 16 types integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Nslip) :: lattice_hex_interactionSlipSlip = & reshape((/& - 1,2,2,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8, & - 2,1,2,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8, & - 2,2,1,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8, & - 6,6,6,1,3,3,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10, & - 6,6,6,3,1,3,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10, & - 6,6,6,3,3,1,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10, & - 7,7,7,9,9,9,1,4,4,4,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & - 7,7,7,9,9,9,4,1,4,4,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & - 7,7,7,9,9,9,4,4,1,4,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & - 7,7,7,9,9,9,4,4,4,1,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & - 7,7,7,9,9,9,4,4,4,4,1,4,11,11,11,11,11,11,11,11,11,11,11,11, & - 7,7,7,9,9,9,4,4,4,4,4,1,11,11,11,11,11,11,11,11,11,11,11,11, & - 8,8,8,10,10,10,11,11,11,11,11,11,1,5,1,5,5,5,5,5,5,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,1,5,1,5,5,5,5,5,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,1,5,1,5,5,5,5,5,5,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,1,5,1,5,5,5,5,5,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,1,5,5,5,1,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,1,5,1,5,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,1,5,5,5,1,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,1,5,1,5,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,1,5,5,5,1,5,5,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,5,5,5,1,5,1, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,1,5,5,5,1,5, & - 8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,5,5,5,1,5,1 & + 1, 5, 5, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14,14,14,14,14,14,14, & + 5, 1, 5, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14,14,14,14,14,14,14, & + 5, 5, 1, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14,14,14,14,14,14,14, & + 15,15,15, 2, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13,13,13,13,13,13,13, & + 15,15,15, 6, 2, 6,10,10,10,10,10,10,13,13,13,13,13,13,13,13,13,13,13,13, & + 15,15,15, 6, 6, 2,10,10,10,10,10,10,13,13,13,13,13,13,13,13,13,13,13,13, & + 18,18,18,16,16,16, 3, 7, 7, 7, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, & + 18,18,18,16,16,16, 7, 3, 7, 7, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, & + 18,18,18,16,16,16, 7, 7, 3, 7, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, & + 18,18,18,16,16,16, 7, 7, 7, 3, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, & + 18,18,18,16,16,16, 7, 7, 7, 7, 3, 7,11,11,11,11,11,11,11,11,11,11,11,11, & + 18,18,18,16,16,16, 7, 7, 7, 7, 7, 3,11,11,11,11,11,11,11,11,11,11,11,11, & + 20,20,20,19,19,19,17,17,17,17,17,17, 4, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 4, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 4, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 4, 8, 4, 8, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 4, 8, 4, 8, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, 8, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 8, 4, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, & + 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 8, 4 & /),(/lattice_hex_Nslip,lattice_hex_Nslip/)) !* isotropic interaction at the moment integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Ntwin) :: lattice_hex_interactionSlipTwin = & reshape((/& - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 & + 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & + 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & + 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & + 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & + 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & + 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, & + 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, & + 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, & + 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16 & /),(/lattice_hex_Nslip,lattice_hex_Ntwin/)) integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Ntwin) :: lattice_hex_interactionTwinTwin = & reshape((/& - 1,2,2,2,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & - 2,1,2,2,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & - 2,2,1,2,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & - 2,2,2,1,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & - 2,2,2,2,1,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & - 2,2,2,2,2,1,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & - 6,6,6,6,6,6,1,3,3,3,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & - 6,6,6,6,6,6,3,1,3,3,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & - 6,6,6,6,6,6,3,3,1,3,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & - 6,6,6,6,6,6,3,3,3,1,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & - 6,6,6,6,6,6,3,3,3,3,1,3,9,9,9,9,9,9,10,10,10,10,10,10, & - 6,6,6,6,6,6,3,3,3,3,3,1,9,9,9,9,9,9,10,10,10,10,10,10, & - 7,7,7,7,7,7,9,9,9,9,9,9,1,4,4,4,4,4,11,11,11,11,11,11, & - 7,7,7,7,7,7,9,9,9,9,9,9,4,1,4,4,4,4,11,11,11,11,11,11, & - 7,7,7,7,7,7,9,9,9,9,9,9,4,4,1,4,4,4,11,11,11,11,11,11, & - 7,7,7,7,7,7,9,9,9,9,9,9,4,4,4,1,4,4,11,11,11,11,11,11, & - 7,7,7,7,7,7,9,9,9,9,9,9,4,4,4,4,1,4,11,11,11,11,11,11, & - 7,7,7,7,7,7,9,9,9,9,9,9,4,4,4,4,4,1,11,11,11,11,11,11, & - 8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,1,5,5,5,5,5, & - 8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,1,5,5,5,5, & - 8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,1,5,5,5, & - 8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,5,1,5,5, & - 8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,5,5,1,5, & - 8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,5,5,5,1 & + 1, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, & + 5, 1, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, & + 5, 5, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, & + 5, 5, 5, 1, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, & + 5, 5, 5, 5, 1, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, & + 5, 5, 5, 5, 5, 1, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, & + 15,15,15,15,15,15, 2, 6, 6, 6, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, & + 15,15,15,15,15,15, 6, 2, 6, 6, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, & + 15,15,15,15,15,15, 6, 6, 2, 6, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, & + 15,15,15,15,15,15, 6, 6, 6, 2, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, & + 15,15,15,15,15,15, 6, 6, 6, 6, 2, 6,10,10,10,10,10,10,13,13,13,13,13,13, & + 15,15,15,15,15,15, 6, 6, 6, 6, 6, 2,10,10,10,10,10,10,13,13,13,13,13,13, & + 18,18,18,18,18,18,16,16,16,16,16,16, 3, 7, 7, 7, 7, 7,11,11,11,11,11,11, & + 18,18,18,18,18,18,16,16,16,16,16,16, 7, 3, 7, 7, 7, 7,11,11,11,11,11,11, & + 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 3, 7, 7, 7,11,11,11,11,11,11, & + 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 7, 3, 7, 7,11,11,11,11,11,11, & + 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 7, 7, 3, 7,11,11,11,11,11,11, & + 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 7, 7, 7, 3,11,11,11,11,11,11, & + 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 4, 8, 8, 8, 8, 8, & + 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 4, 8, 8, 8, 8, & + 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 4, 8, 8, 8, & + 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 8, 4, 8, 8, & + 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 4, 8, & + 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 4 & /),(/lattice_hex_Ntwin,lattice_hex_Ntwin/)) !* isotropic interaction at the moment integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Nslip) :: lattice_hex_interactionTwinSlip = & reshape((/& - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & - 1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 & + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, & + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, & + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, & + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, & + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, & + 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, & + 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, & + 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, & + 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, & + 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, & + 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, & + 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, & + 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, & + 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, & + 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, & + 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, & + 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, & + 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, & + 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, & + 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, & + 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, & + 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, & + 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, & + 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16 & /),(/lattice_hex_Ntwin,lattice_hex_Nslip/)) @@ -585,42 +593,46 @@ subroutine lattice_init() !************************************** !* Module initialization * !************************************** - use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error - use material, only: material_configfile,material_partPhase - implicit none + use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error + use material, only: material_configfile,material_partPhase + implicit none - integer(pInt), parameter :: fileunit = 200 - integer(pInt) i,Nsections + integer(pInt), parameter :: fileunit = 200 + integer(pInt) i,Nsections - write(6,*) - write(6,*) '<<<+- lattice init -+>>>' - write(6,*) - - if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file - Nsections = IO_countSections(fileunit,material_partPhase) - lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex - close(fileunit) + write(6,*) + write(6,*) '<<<+- lattice init -+>>>' + write(6,*) - allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal - allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal - allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal - allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal - allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal + if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file + Nsections = IO_countSections(fileunit,material_partPhase) + lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex + close(fileunit) - allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal - allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal - allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal - allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal - allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal - allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal - allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal + allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal + allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal + allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal + allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal + allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal - allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt - allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt - allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt - allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt - -endsubroutine + allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal + allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal + allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal + allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal + allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal + allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal + + allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal + + allocate(lattice_NslipSystem(lattice_maxNslipFamily,lattice_Nstructure)); lattice_NslipSystem = 0.0_pReal + allocate(lattice_NtwinSystem(lattice_maxNtwinFamily,lattice_Nstructure)); lattice_NtwinSystem = 0.0_pReal + + allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt + allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt + allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt + allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt + +end subroutine function lattice_initializeStructure(struct,CoverA) @@ -643,14 +655,20 @@ function lattice_initializeStructure(struct,CoverA) real(pReal), dimension(lattice_maxNtwin) :: ts = 0.0_pReal real(pReal), dimension(3) :: hex_d = 0.0_pReal, & hex_n = 0.0_pReal + integer(pInt), dimension(lattice_maxNslipFamily) :: myNslipSystem = 0_pInt + integer(pInt), dimension(lattice_maxNtwinFamily) :: myNtwinSystem = 0_pInt integer(pInt) :: i,myNslip,myNtwin,myStructure = 0_pInt - logical :: processMe = .false. + logical :: processMe integer(pInt) lattice_initializeStructure + processMe = .false. + select case(struct(1:3)) ! check first three chars of structure name case ('fcc') myStructure = 1_pInt + myNslipSystem = lattice_fcc_NslipSystem + myNtwinSystem = lattice_fcc_NtwinSystem myNslip = lattice_fcc_Nslip myNtwin = lattice_fcc_Ntwin lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt @@ -675,6 +693,8 @@ function lattice_initializeStructure(struct,CoverA) case ('bcc') myStructure = 2_pInt + myNslipSystem = lattice_bcc_NslipSystem + myNtwinSystem = lattice_bcc_NtwinSystem myNslip = lattice_bcc_Nslip myNtwin = lattice_bcc_Ntwin lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt @@ -701,29 +721,31 @@ function lattice_initializeStructure(struct,CoverA) if (CoverA > 0.0_pReal) then lattice_hex_Nstructure = lattice_hex_Nstructure + 1_pInt myStructure = 2_pInt + lattice_hex_Nstructure + myNslipSystem = lattice_hex_NslipSystem + myNtwinSystem = lattice_hex_NtwinSystem myNslip = lattice_hex_Nslip myNtwin = lattice_hex_Ntwin processMe = .true. ! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c) do i = 1,myNslip - hex_d(1) = lattice_hex_systemSlip(1,i) ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] - hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))/dsqrt(3.0_pReal) - hex_d(3) = lattice_hex_systemSlip(4,i)/CoverA - hex_n(1) = lattice_hex_systemSlip(5,i)*1.5_pReal ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) - hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))*(0.5_pReal*dsqrt(3.0_pReal)) - hex_n(3) = lattice_hex_systemSlip(8,i)*CoverA + hex_d(1) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] + hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*dsqrt(3.0_pReal)) + hex_d(3) = lattice_hex_systemSlip(4,i)*CoverA + hex_n(1) = lattice_hex_systemSlip(5,i) ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) + hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/dsqrt(3.0_pReal) + hex_n(3) = lattice_hex_systemSlip(8,i)/CoverA sd(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d)) sn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n)) st(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) enddo do i = 1,myNtwin - hex_d(1) = lattice_hex_systemTwin(1,i) - hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))/dsqrt(3.0_pReal) - hex_d(3) = lattice_hex_systemTwin(4,i)/CoverA - hex_n(1) = lattice_hex_systemTwin(5,i)*1.5_pReal - hex_n(2) = (lattice_hex_systemTwin(5,i)+2.0_pReal*lattice_hex_systemTwin(6,i))*(0.5_pReal*dsqrt(3.0_pReal)) - hex_n(3) = lattice_hex_systemTwin(8,i)*CoverA + hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal + hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*dsqrt(3.0_pReal)) + hex_d(3) = lattice_hex_systemTwin(4,i)*CoverA + hex_n(1) = lattice_hex_systemTwin(4,i) + hex_n(2) = (lattice_hex_systemTwin(4,i)+2.0_pReal*lattice_hex_systemTwin(6,i))/dsqrt(3.0_pReal) + hex_n(3) = lattice_hex_systemTwin(8,i)/CoverA td(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d)) tn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n)) @@ -754,15 +776,17 @@ function lattice_initializeStructure(struct,CoverA) lattice_Qtwin(:,:,i,myStructure) = math_RodrigToR(tn(:,i),180.0_pReal*inRad) lattice_shearTwin(i,myStructure) = ts(i) enddo + lattice_NslipSystem(1:lattice_maxNslipFamily,myStructure) = myNslipSystem ! number of slip systems in each family + lattice_NtwinSystem(1:lattice_maxNtwinFamily,myStructure) = myNtwinSystem ! number of twin systems in each family lattice_interactionSlipSlip(1:myNslip,1:myNslip,myStructure) = interactionSlipSlip(1:myNslip,1:myNslip) lattice_interactionSlipTwin(1:myNslip,1:myNtwin,myStructure) = interactionSlipTwin(1:myNslip,1:myNtwin) - lattice_interactionTwinTwin(1:myNtwin,1:myNtwin,myStructure) = interactionTwinTwin(1:myNtwin,1:myNtwin) lattice_interactionTwinSlip(1:myNtwin,1:myNslip,myStructure) = interactionTwinSlip(1:myNtwin,1:myNslip) + lattice_interactionTwinTwin(1:myNtwin,1:myNtwin,myStructure) = interactionTwinTwin(1:myNtwin,1:myNtwin) endif lattice_initializeStructure = myStructure -endfunction +end function END MODULE diff --git a/code/material.config b/code/material.config index ec7bf1d7b..083fcf441 100644 --- a/code/material.config +++ b/code/material.config @@ -61,58 +61,70 @@ Ngrains 100 ##################### -[Aluminum] # below given format will not work. need to select one constitution block from it. +[Aluminum_J2isotropic] + +constitution j2 + +(output) flowstress +(output) strainrate + +c11 110.9e9 +c12 58.34e9 +taylorfactor 3 +tau0 31e6 +gdot0 0.001 +n 20 +h0 75e6 +tausat 63e6 +w0 1 + +[Aluminum_phenopowerlaw] +# slip only +constitution phenopowerlaw + +(output) resistance_slip +(output) shearrate_slip +(output) resolvedstress_slip +(output) totalshear +(output) resistance_twin +(output) shearrate_twin +(output) resolvedstress_twin +(output) totalvolfrac + +lattice_structure fcc +Nslip 12 0 0 0 # per family +Ntwin 0 0 0 0 # per family c11 106.75e9 c12 60.41e9 c44 28.34e9 -lattice_structure fcc -Nslip 12 -Ntwin 12 -constitution phenomenological -(output) slipResistance -(output) rateOfShear -### phenomenological constitutive parameters ### -s0_slip 31e6 gdot0_slip 0.001 n_slip 20 -h0 75e6 -s_sat 63e6 -w0 2.25 -latent_ratio 1.4 +s0_slip 31e6 0 0 0 # per family +ssat_slip 63e6 0 0 0 # per family +gdot0_twin 0.001 +n_twin 20 +s0_twin 31e6 0 0 0 # per family +s_pr 0 # push-up factor for slip saturation due to twinning +twin_b 0 +twin_c 0 +twin_d 0 +twin_e 0 +h0_slipslip 75e6 +h0_sliptwin 0 +h0_twinslip 0 +h0_twintwin 0 +interaction_slipslip 1 1 1.4 1.4 1.4 1.4 +interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 +interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 +interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -#constitution j2 -#c11 110.9e9 -#c12 58.34e9 -#(output) flowstress -#(output) strainrate -#taylorfactor 3 -#s0 31e6 -#gdot0 0.001 -#n 20 -#h0 75e6 -#s_sat 63e6 -#w0 2.25 - -[Copper] - -[Ferrite] - -[Martensite] [TWIP steel FeMnC] -C11 175.0e9 # elastic constants in Pa -C12 115.0e9 -C44 135.0e9 -lattice_structure fcc -Nslip 12 -Ntwin 12 -#constitution phenomenological -#(output) slipResistance -#(output) rateOfShear constitution dislobased + (output) state_slip (output) rateofshear_slip (output) mfp_slip @@ -122,14 +134,12 @@ constitution dislobased (output) mfp_twin (output) thresholdstress_twin -### phenomenological constitutive parameters ### -#s0_slip 85.0e6 # initial slip resistance -#gdot0_slip 0.001 # reference shear rate -#n_slip 100.0 # stress exponent -#h0 355.0e6 # initial hardening slope -#s_sat 265.0e6 # saturation stress -#w0 1.0 # exponent -#latent_ratio 1.4 # latent/self hardening ratio +C11 175.0e9 # elastic constants in Pa +C12 115.0e9 +C44 135.0e9 +lattice_structure fcc +Nslip 12 +Ntwin 12 ### dislocation density-based constitutive parameters ### burgers 2.56e-10 # Burgers vector [m] diff --git a/code/material.f90 b/code/material.f90 index 8d6f7bd92..8da46dc94 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -28,7 +28,7 @@ integer(pInt) material_Nhomogenization, & material_Nphase, & ! number of phases material_Ntexture, & ! number of textures microstructure_maxNconstituents, & ! max number of constituents in any phase - homogenization_maxNgrains, & ! max number of grains in any homogenization + homogenization_maxNgrains, & ! max number of grains in any USED homogenization texture_maxNgauss, & ! max number of Gauss components in any texture texture_maxNfiber ! max number of Fiber components in any texture character(len=64), dimension(:), allocatable :: homogenization_name, & ! name of each homogenization @@ -44,10 +44,12 @@ integer(pInt), dimension(:), allocatable :: homogenization_Ngrains, & microstructure_Nconstituents, & ! number of constituents in each microstructure phase_constitutionInstance, & ! instance of particular constitution of each phase phase_Noutput, & ! number of '(output)' items per phase - phase_localConstitution, & ! flag phases with local constitutive law texture_symmetry, & ! number of symmetric orientations per texture texture_Ngauss, & ! number of Gauss components per texture texture_Nfiber ! number of Fiber components per texture +logical, dimension(:), allocatable :: homogenization_active, & ! + microstructure_active, & ! + phase_localConstitution ! flags phases with local constitutive law integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase IDs of each microstructure microstructure_texture ! texture IDs of each microstructure real(pReal), dimension(:,:), allocatable :: microstructure_fraction ! vol fraction of each constituent in microstructure @@ -71,7 +73,7 @@ subroutine material_init() !* Definition of variables integer(pInt), parameter :: fileunit = 200 - integer(pInt) i + integer(pInt) i,j write(6,*) write(6,*) '<<<+- material init -+>>>' @@ -94,14 +96,22 @@ subroutine material_init() write (6,*) write (6,*) 'MATERIAL configuration' write (6,*) - write (6,*) 'Homogenization' + write (6,'(a32,x,a16,x,a6)') 'homogenization ','type ','grains' do i = 1,material_Nhomogenization - write (6,'(a32,x,a16,x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i) + write (6,'(x,a32,x,a16,x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i) enddo write (6,*) - write (6,*) 'Microstructure' + write (6,'(a32,x,a12)') 'microstructure ','constituents' do i = 1,material_Nmicrostructure write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i) + if (microstructure_Nconstituents(i) > 0_pInt) then + do j = 1,microstructure_Nconstituents(i) + write (6,'(a1,x,a32,x,a32,x,f6.4)') '>',phase_name(microstructure_phase(j,i)),& + texture_name(microstructure_texture(j,i)),& + microstructure_fraction(j,i) + enddo + write (6,*) + endif enddo call material_populateGrains() @@ -115,6 +125,7 @@ subroutine material_parseHomogenization(file,myPart) use prec, only: pInt use IO + use mesh, only: mesh_element implicit none character(len=*), intent(in) :: myPart @@ -128,16 +139,15 @@ subroutine material_parseHomogenization(file,myPart) Nsections = IO_countSections(file,myPart) material_Nhomogenization = Nsections - write (6,*) 'homogenization sections found',material_Nhomogenization allocate(homogenization_name(Nsections)); homogenization_name = '' allocate(homogenization_type(Nsections)); homogenization_type = '' allocate(homogenization_typeInstance(Nsections)); homogenization_typeInstance = 0_pInt allocate(homogenization_Ngrains(Nsections)); homogenization_Ngrains = 0_pInt allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt - - write(6,*) 'scanning for (output)',homogenization_Noutput + allocate(homogenization_active(Nsections)); homogenization_active = .false. + + forall (s = 1:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) - write(6,*) 'count of (output)',homogenization_Noutput rewind(file) line = '' @@ -171,7 +181,7 @@ subroutine material_parseHomogenization(file,myPart) endif enddo -100 homogenization_maxNgrains = maxval(homogenization_Ngrains) +100 homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active) return endsubroutine @@ -183,6 +193,7 @@ subroutine material_parseMicrostructure(file,myPart) use prec, only: pInt use IO + use mesh, only: mesh_element implicit none character(len=*), intent(in) :: myPart @@ -197,6 +208,9 @@ subroutine material_parseMicrostructure(file,myPart) material_Nmicrostructure = Nsections allocate(microstructure_name(Nsections)); microstructure_name = '' allocate(microstructure_Nconstituents(Nsections)) + allocate(microstructure_active(Nsections)) + + forall (i = 1:Nsections) microstructure_active(i) = any(mesh_element(4,:) == i) ! current microstructure used in model? microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections) microstructure_maxNconstituents = maxval(microstructure_Nconstituents) @@ -452,7 +466,7 @@ subroutine material_populateGrains() allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure)); Ngrains = 0_pInt -! count grains per homog/micro pair +! identify maximum grain count per IP (from element) and find grains per homog/micro pair do e = 1,mesh_NcpElems homog = mesh_element(3,e) micro = mesh_element(4,e) @@ -463,7 +477,7 @@ subroutine material_populateGrains() Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e)) enddo - allocate(volumeOfGrain(maxval(Ngrains))) ! reserve memory for maximum case + allocate(volumeOfGrain(maxval(Ngrains))) ! reserve memory for maximum case allocate(phaseOfGrain(maxval(Ngrains))) ! reserve memory for maximum case allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case diff --git a/code/mpie_cpfem_marc.f90 b/code/mpie_cpfem_marc.f90 index ac7f9b353..a6e4ef42d 100644 --- a/code/mpie_cpfem_marc.f90 +++ b/code/mpie_cpfem_marc.f90 @@ -43,7 +43,7 @@ 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_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug + include "constitutive_phenopowerlaw.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.f90" ! uses prec, IO, math, lattice, mesh, debug diff --git a/code/numerics.f90 b/code/numerics.f90 index f40aea921..cb04b87a9 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -16,8 +16,8 @@ integer(pInt) iJacoStiffness, & ! freque real(pReal) relevantStrain, & ! strain increment considered significant pert_Fg, & ! strain perturbation for FEM Jacobi subStepMin, & ! minimum (relative) size of sub-step allowed during cutback in crystallite - rTol_crystalliteState, & ! relative tolerance in crystallite state loop - rTol_crystalliteTemperature, & ! relative tolerance in crystallite temperature loop + rTol_crystalliteState, & ! relative tolerance in crystallite state loop + rTol_crystalliteTemperature, & ! relative tolerance in crystallite temperature loop rTol_crystalliteStress, & ! relative tolerance in crystallite stress loop aTol_crystalliteStress, & ! absolute tolerance in crystallite stress loop resToler, & ! relative tolerance of residual in GIA iteration @@ -57,25 +57,6 @@ subroutine numerics_init() character(len=64) tag character(len=1024) line - !*** global variables ***! - ! relevantStrain - ! iJacoStiffness - ! iJacoLpresiduum - ! pert_Fg - ! nHomog - ! nCryst - ! nState - ! nStress - ! subStepMin - ! rTol_crystalliteState - ! rTol_crystalliteTemperature - ! rTol_crystalliteStress - ! aTol_crystalliteStress - ! resToler - ! resAbsol - ! resBound - ! NRiterMax - write(6,*) write(6,*) '<<<+- numerics init -+>>>' write(6,*) @@ -90,8 +71,8 @@ subroutine numerics_init() nState = 10_pInt nStress = 40_pInt subStepMin = 1.0e-3_pReal - rTol_crystalliteState = 1.0e-6_pReal - rTol_crystalliteTemperature = 1.0e-6_pReal + rTol_crystalliteState = 1.0e-6_pReal + rTol_crystalliteTemperature = 1.0e-6_pReal rTol_crystalliteStress = 1.0e-6_pReal aTol_crystalliteStress = 1.0e-8_pReal resToler = 1.0e-4_pReal @@ -132,8 +113,8 @@ subroutine numerics_init() case ('substepmin') subStepMin = IO_floatValue(line,positions,2) case ('rtol_crystallitestate') - rTol_crystalliteState = IO_floatValue(line,positions,2) - case ('rtol_crystallitetemperature') + rTol_crystalliteState = IO_floatValue(line,positions,2) + case ('rtol_crystallitetemperature') rTol_crystalliteTemperature = IO_floatValue(line,positions,2) case ('rtol_crystallitestress') rTol_crystalliteStress = IO_floatValue(line,positions,2) @@ -169,7 +150,7 @@ subroutine numerics_init() write(6,'(a24,x,i8)') 'nState: ',nState write(6,'(a24,x,i8)') 'nStress: ',nStress write(6,'(a24,x,e8.1)') 'subStepMin: ',subStepMin - write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState + write(6,'(a24,x,e8.1)') 'rTol_crystalliteState: ',rTol_crystalliteState write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress