From d6abc002187c29aafe7cde49d324c937ec2f9a36 Mon Sep 17 00:00:00 2001 From: Chen Zhang Date: Tue, 13 Oct 2015 18:52:01 +0000 Subject: [PATCH] add pheno+ module in --- code/Makefile | 137 ++-- code/constitutive.f90 | 222 +++--- code/crystallite.f90 | 262 +++--- code/debug.f90 | 60 +- code/material.f90 | 243 +++--- code/plastic_phenoplus.f90 | 1356 ++++++++++++++++++++++++++++++++ code/plastic_phenopowerlaw.f90 | 102 +-- 7 files changed, 1899 insertions(+), 483 deletions(-) create mode 100644 code/plastic_phenoplus.f90 diff --git a/code/Makefile b/code/Makefile index 83dd3ff89..a088f0bad 100644 --- a/code/Makefile +++ b/code/Makefile @@ -93,7 +93,7 @@ OPTIMIZATION_OFF_gfortran :=-O0 OPTIMIZATION_DEFENSIVE_ifort :=-O2 OPTIMIZATION_DEFENSIVE_gfortran :=-O2 OPTIMIZATION_AGGRESSIVE_ifort :=-O3 $(PORTABLE_SWITCH) -no-prec-div -fp-model fast=2 -ipo -OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize +OPTIMIZATION_AGGRESSIVE_gfortran :=-O3 $(PORTABLE_SWITCH) -ffast-math -funroll-loops -ftree-vectorize LINK_OPTIONS_ifort :=-shared-intel @@ -119,7 +119,7 @@ endif #-ftz: flush unterflow to zero, automatically set if O<0,1,2,3> >0 #-assume byterecl record length is given in bytes (also set by -standard-semantics) #-fimplicit-none: assume "implicit-none" even if not present in source -#-diag-disable: disables warnings, where +#-diag-disable: disables warnings, where # warning ID 5268: the text exceeds right hand column allowed on the line (we have only comments there) #-warn: enables warnings, where # declarations: any undeclared names (alternative name: -implicitnone) @@ -154,10 +154,10 @@ DEBUG_OPTIONS_ifort :=-g\ #COMPILE SWITCHES FOR RUNTIME DEBUGGING #-g: Generate symbolic debugging information in the object file #-traceback: Generate extra information in the object file to provide source file traceback information when a severe error occurs at run time. -#-gen-interfaces: Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/ +#-gen-interfaces: Generate an interface block for each routine. http://software.intel.com/en-us/blogs/2012/01/05/doctor-fortran-gets-explicit-again/ #-fp-stack-check: Generate extra code after every function call to ensure that the floating-point (FP) stack is in the expected state. #-ftrapuv Trap uninitalized variables -#-check: checks at runtime, where +#-check: checks at runtime, where # bounds: check if an array index is too small (<1) or too large! # format: Checking for the data type of an item being formatted for output. # output_conversion: Checking for the fit of data items within a designated format descriptor field. @@ -171,7 +171,7 @@ DEBUG_OPTIONS_ifort :=-g\ ################################################################################################### #MORE OPTIONS FOR RUNTIME DEBUGGING #-heap-arrays: should not be done for OpenMP, but set "ulimit -s unlimited" on shell. Probably it helps also to unlimit other limits -#-check: checks at runtime, where +#-check: checks at runtime, where # arg_temp_created: will cause a lot of warnings because we create a bunch of temporary arrays (performance?) # stack: LINK_OPTIONS_gfortran :=-Wl,-undefined,dynamic_lookup @@ -220,48 +220,48 @@ endif # -Wreal-q-constant: warn about real-literal-constants with 'q' exponent-letter # -Wunused: a number of unused-xxx warnings # these are general (non -Fortran options) implied by -Wall -# -Waddress -# -Warray-bounds (only with -O2) -# -Wc++11-compat -# -Wchar-subscripts -# -Wcomment -# -Wformat -# -Wmaybe-uninitialized -# -Wnonnull -# -Wparentheses -# -Wpointer-sign -# -Wreorder -# -Wreturn-type -# -Wsequence-point -# -Wstrict-aliasing -# -Wstrict-overflow=1 -# -Wswitch -# -Wtrigraphs -# -Wuninitialized -# -Wunknown-pragmas -# -Wunused-function -# -Wunused-label -# -Wunused-value -# -Wunused-variable -# -Wvolatile-register-var +# -Waddress +# -Warray-bounds (only with -O2) +# -Wc++11-compat +# -Wchar-subscripts +# -Wcomment +# -Wformat +# -Wmaybe-uninitialized +# -Wnonnull +# -Wparentheses +# -Wpointer-sign +# -Wreorder +# -Wreturn-type +# -Wsequence-point +# -Wstrict-aliasing +# -Wstrict-overflow=1 +# -Wswitch +# -Wtrigraphs +# -Wuninitialized +# -Wunknown-pragmas +# -Wunused-function +# -Wunused-label +# -Wunused-value +# -Wunused-variable +# -Wvolatile-register-var #-Wextra: sets the following Fortran options: # -Wunuses-parameter: # -Wcompare-reals: # these are general (non -Fortran options) implied by -Wextra -# -Wclobbered -# -Wempty-body -# -Wignored-qualifiers -# -Wmissing-field-initializers -# -Woverride-init -# -Wsign-compare -# -Wtype-limits -# -Wuninitialized +# -Wclobbered +# -Wempty-body +# -Wignored-qualifiers +# -Wmissing-field-initializers +# -Woverride-init +# -Wsign-compare +# -Wtype-limits +# -Wuninitialized # -Wunused-but-set-parameter (only with -Wunused or -Wall) # -Wno-globals ################################################################################################### #MORE OPTIONS FOR DEBUGGING DURING COMPILATION -#-Warray-temporarieswarnings: because we have many temporary arrays (performance issue?): +#-Warray-temporarieswarnings: because we have many temporary arrays (performance issue?): #-Wimplicit-interface: no interfaces for lapack routines #-Wunsafe-loop-optimizations: warn if the loop cannot be optimized due to nontrivial assumptions. #-Wstrict-overflow: @@ -310,16 +310,16 @@ SOURCE_FILES = \ source_thermal_dissipation.o source_thermal_externalheat.o \ source_damage_isoBrittle.o source_damage_isoDuctile.o source_damage_anisoBrittle.o source_damage_anisoDuctile.o \ source_vacancy_phenoplasticity.o source_vacancy_irradiation.o source_vacancy_thermalfluc.o - + KINEMATICS_FILES = \ kinematics_cleavage_opening.o kinematics_slipplane_opening.o \ kinematics_thermal_expansion.o \ - kinematics_vacancy_strain.o kinematics_hydrogen_strain.o - + kinematics_vacancy_strain.o kinematics_hydrogen_strain.o + PLASTIC_FILES = \ plastic_dislotwin.o plastic_disloUCLA.o plastic_disloKMC.o plastic_j2.o plastic_phenopowerlaw.o \ - plastic_titanmod.o plastic_nonlocal.o plastic_none.o - + plastic_titanmod.o plastic_nonlocal.o plastic_none.o plastic_phenoplus.o + THERMAL_FILES = \ thermal_isothermal.o thermal_adiabatic.o thermal_conduction.o @@ -336,7 +336,7 @@ HYDROGENFLUX_FILES = \ hydrogenflux_isoconc.o hydrogenflux_cahnhilliard.o HOMOGENIZATION_FILES = \ - homogenization_RGC.o homogenization_isostrain.o homogenization_none.o + homogenization_RGC.o homogenization_isostrain.o homogenization_none.o ##################### # Spectral Solver @@ -367,31 +367,31 @@ SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o DAMASK_spectral_utilities.o \ $(SPECTRAL_SOLVER_FILES) -DAMASK_spectral.exe: DAMASK_spectral_driver.o +DAMASK_spectral.exe: DAMASK_spectral_driver.o $(PREFIX) $(LINKERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \ -o DAMASK_spectral.exe DAMASK_spectral_driver.o \ $(SPECTRAL_FILES) $(LIBRARIES) $(SUFFIX) - - + + DAMASK_spectral_driver.o: DAMASK_spectral_driver.f90 \ $(SPECTRAL_SOLVER_FILES) $(PREFIX) $(COMPILERNAME) $(COMPILE_MAXOPTI) -c DAMASK_spectral_driver.f90 $(SUFFIX) DAMASK_spectral_solverAL.o: DAMASK_spectral_solverAL.f90 \ DAMASK_spectral_utilities.o - + DAMASK_spectral_solverPolarisation.o: DAMASK_spectral_solverPolarisation.f90 \ DAMASK_spectral_utilities.o - + DAMASK_spectral_solverBasicPETSc.o: DAMASK_spectral_solverBasicPETSc.f90 \ DAMASK_spectral_utilities.o spectral_thermal.o: spectral_thermal.f90 \ DAMASK_spectral_utilities.o - + spectral_damage.o: spectral_damage.f90 \ DAMASK_spectral_utilities.o - + DAMASK_spectral_utilities.o: DAMASK_spectral_utilities.f90 \ CPFEM.o @@ -443,7 +443,7 @@ FEM_hydrogenflux.o: FEM_hydrogenflux.f90 \ FEM_utilities.o FEM_utilities.o: FEM_utilities.f90 \ - CPFEM.o + CPFEM.o FEZoo.o: $(wildcard FEZoo.f90) \ IO.o @@ -513,8 +513,8 @@ homogenization_none.o: homogenization_none.f90 \ crystallite.o crystallite.o: crystallite.f90 \ - constitutive.o - + constitutive.o + constitutive.o: constitutive.f90 \ $(SOURCE_FILES) \ $(KINEMATICS_FILES) \ @@ -531,13 +531,13 @@ source_damage_isoBrittle.o: source_damage_isoBrittle.f90 \ source_damage_isoDuctile.o: source_damage_isoDuctile.f90 \ lattice.o - + source_damage_anisoBrittle.o: source_damage_anisoBrittle.f90 \ lattice.o - + source_damage_anisoDuctile.o: source_damage_anisoDuctile.f90 \ lattice.o - + source_vacancy_phenoplasticity.o: source_vacancy_phenoplasticity.f90 \ lattice.o @@ -546,22 +546,22 @@ source_vacancy_irradiation.o: source_vacancy_irradiation.f90 \ source_vacancy_thermalfluc.o: source_vacancy_thermalfluc.f90 \ lattice.o - + kinematics_cleavage_opening.o: kinematics_cleavage_opening.f90 \ lattice.o - + kinematics_slipplane_opening.o: kinematics_slipplane_opening.f90 \ lattice.o - + kinematics_thermal_expansion.o: kinematics_thermal_expansion.f90 \ lattice.o - + kinematics_vacancy_strain.o: kinematics_vacancy_strain.f90 \ lattice.o - + kinematics_hydrogen_strain.o: kinematics_hydrogen_strain.f90 \ lattice.o - + plastic_nonlocal.o: plastic_nonlocal.f90 \ lattice.o @@ -580,6 +580,9 @@ plastic_dislotwin.o: plastic_dislotwin.f90 \ plastic_phenopowerlaw.o: plastic_phenopowerlaw.f90 \ lattice.o +plastic_phenoplus.o: plastic_phenoplus.f90 \ + lattice.o + plastic_j2.o: plastic_j2.f90 \ lattice.o @@ -616,7 +619,7 @@ debug.o: debug.f90 \ numerics.o: numerics.f90 \ libs.o - + libs.o: libs.f90 \ IO.o @@ -628,17 +631,17 @@ DAMASK_interface.o: DAMASK_spectral_interface.f90 \ $(wildcard DAMASK_FEM_interface.f90) \ prec.o $(PREFIX) $(COMPILERNAME) $(COMPILE) -c $(INTERFACENAME) -fall-intrinsics -o DAMASK_interface.o $(SUFFIX) -#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored +#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored # and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external # --> allows the use of 'getcwd' prec.o: prec.f90 $(PREFIX) $(COMPILERNAME) $(COMPILE) -c prec.f90 -fno-range-check -fall-intrinsics -fno-fast-math $(SUFFIX) # fno-range-check: Disable range checking on results of simplification of constant expressions during compilation # --> allows the definition of DAMASK_NaN -#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored +#-fall-intrinsics: all intrinsic procedures (including the GNU-specific extensions) are accepted. -Wintrinsics-std will be ignored # and no user-defined procedure with the same name as any intrinsic will be called except when it is explicitly declared external # --> allows the use of 'isnan' -#-fno-fast-math: +#-fno-fast-math: # --> otherwise, when setting -ffast-math, isnan always evaluates to false (I would call it a bug) else DAMASK_interface.o: DAMASK_spectral_interface.f90 \ diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 4a17c4b36..5c7afd5c4 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -8,7 +8,7 @@ module constitutive use prec, only: & pInt - + implicit none private integer(pInt), public, protected :: & @@ -17,7 +17,7 @@ module constitutive constitutive_source_maxSizePostResults, & constitutive_source_maxSizeDotState - public :: & + public :: & constitutive_init, & constitutive_homogenizedC, & constitutive_microstructure, & @@ -28,10 +28,10 @@ module constitutive constitutive_collectDotState, & constitutive_collectDeltaState, & constitutive_postResults - + private :: & constitutive_hooke_TandItsTangent - + contains @@ -66,8 +66,8 @@ subroutine constitutive_init() use material, only: & material_phase, & material_Nphase, & - material_localFileExt, & - material_configFile, & + material_localFileExt, & + material_configFile, & phase_name, & phase_plasticity, & phase_plasticityInstance, & @@ -78,6 +78,7 @@ subroutine constitutive_init() PLASTICITY_none_ID, & PLASTICITY_j2_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & PLASTICITY_disloucla_ID, & @@ -101,6 +102,7 @@ subroutine constitutive_init() PLASTICITY_NONE_label, & PLASTICITY_J2_label, & PLASTICITY_PHENOPOWERLAW_label, & + PLASTICITY_PHENOPLUS_label, & PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOKMC_label, & PLASTICITY_DISLOUCLA_label, & @@ -116,11 +118,12 @@ subroutine constitutive_init() SOURCE_vacancy_irradiation_label, & SOURCE_vacancy_thermalfluc_label, & plasticState, & - sourceState + sourceState use plastic_none use plastic_j2 use plastic_phenopowerlaw + use plastic_phenoplus use plastic_dislotwin use plastic_dislokmc use plastic_disloucla @@ -155,7 +158,7 @@ subroutine constitutive_init() character(len=32) :: outputName !< name of output, intermediate fix until HDF5 output is ready logical :: knownPlasticity, knownSource, nonlocalConstitutionPresent nonlocalConstitutionPresent = .false. - + !-------------------------------------------------------------------------------------------------- ! parse plasticities from config file if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... @@ -163,6 +166,7 @@ subroutine constitutive_init() if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_J2_ID)) call plastic_j2_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT) + if (any(phase_plasticity == PLASTICITY_PHENOPLUS_ID)) call plastic_phenoplus_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOKMC_ID)) call plastic_dislokmc_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init(FILEUNIT) @@ -187,7 +191,7 @@ subroutine constitutive_init() if (any(phase_source == SOURCE_vacancy_irradiation_ID)) call source_vacancy_irradiation_init(FILEUNIT) if (any(phase_source == SOURCE_vacancy_thermalfluc_ID)) call source_vacancy_thermalfluc_init(FILEUNIT) close(FILEUNIT) - + !-------------------------------------------------------------------------------------------------- ! parse kinematic mechanisms from config file if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present... @@ -199,17 +203,17 @@ subroutine constitutive_init() if (any(phase_kinematics == KINEMATICS_hydrogen_strain_ID)) call kinematics_hydrogen_strain_init(FILEUNIT) close(FILEUNIT) - mainProcess: if (worldrank == 0) then + mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- constitutive init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" endif mainProcess - + !-------------------------------------------------------------------------------------------------- ! write description file for constitutive phase output if (worldrank == 0_pInt) then - call IO_write_jobFile(FILEUNIT,'outputConstitutive') + call IO_write_jobFile(FILEUNIT,'outputConstitutive') do phase = 1_pInt,material_Nphase if (any(material_phase == phase)) then ! is this phase active? instance = phase_plasticityInstance(phase) ! which instance of a plasticity is present phase @@ -230,6 +234,11 @@ subroutine constitutive_init() thisNoutput => plastic_phenopowerlaw_Noutput thisOutput => plastic_phenopowerlaw_output thisSize => plastic_phenopowerlaw_sizePostResult + case (PLASTICITY_PHENOPLUS_ID) + outputName = PLASTICITY_PHENOPLUS_label + thisNoutput => plastic_phenoplus_Noutput + thisOutput => plastic_phenoplus_output + thisSize => plastic_phenoplus_sizePostResult case (PLASTICITY_DISLOTWIN_ID) outputName = PLASTICITY_DISLOTWIN_label thisNoutput => plastic_dislotwin_Noutput @@ -244,7 +253,7 @@ subroutine constitutive_init() outputName = PLASTICITY_DISLOUCLA_label thisNoutput => plastic_disloucla_Noutput thisOutput => plastic_disloucla_output - thisSize => plastic_disloucla_sizePostResult + thisSize => plastic_disloucla_sizePostResult case (PLASTICITY_TITANMOD_ID) outputName = PLASTICITY_TITANMOD_label thisNoutput => plastic_titanmod_Noutput @@ -257,7 +266,7 @@ subroutine constitutive_init() thisSize => plastic_nonlocal_sizePostResult case default knownPlasticity = .false. - end select + end select write(FILEUNIT,'(/,a,/)') '['//trim(phase_name(phase))//']' if (knownPlasticity) then write(FILEUNIT,'(a)') '(plasticity)'//char(9)//trim(outputName) @@ -324,21 +333,21 @@ subroutine constitutive_init() thisNoutput => source_vacancy_thermalfluc_Noutput thisOutput => source_vacancy_thermalfluc_output thisSize => source_vacancy_thermalfluc_sizePostResult - case default + case default knownSource = .false. - end select + end select if (knownSource) then write(FILEUNIT,'(a)') '(source)'//char(9)//trim(outputName) do e = 1_pInt,thisNoutput(instance) write(FILEUNIT,'(a,i4)') trim(thisOutput(e,instance))//char(9),thisSize(e,instance) enddo endif - enddo + enddo endif enddo close(FILEUNIT) endif - + constitutive_plasticity_maxSizeDotState = 0_pInt constitutive_plasticity_maxSizePostResults = 0_pInt constitutive_source_maxSizeDotState = 0_pInt @@ -351,7 +360,7 @@ subroutine constitutive_init() sourceState(phase)%p(mySource)%partionedState0 = sourceState(phase)%p(mySource)%State0 forall(mySource = 1_pInt:phase_Nsources(phase)) & sourceState(phase)%p(mySource)%State = sourceState(phase)%p(mySource)%State0 - + constitutive_plasticity_maxSizeDotState = max(constitutive_plasticity_maxSizeDotState, & plasticState(phase)%sizeDotState) constitutive_plasticity_maxSizePostResults = max(constitutive_plasticity_maxSizePostResults, & @@ -378,7 +387,7 @@ subroutine constitutive_init() ! report constitutive_maxSizeState = maxval(constitutive_sizeState) constitutive_plasticity_maxSizeDotState = maxval(constitutive_sizeDotState) - + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then write(6,'(a32,1x,7(i8,1x))') 'constitutive_state0: ', shape(constitutive_state0) write(6,'(a32,1x,7(i8,1x))') 'constitutive_partionedState0: ', shape(constitutive_partionedState0) @@ -406,7 +415,7 @@ end subroutine constitutive_init !-------------------------------------------------------------------------------------------------- function constitutive_homogenizedC(ipc,ip,el) use prec, only: & - pReal + pReal use material, only: & phase_plasticity, & material_phase, & @@ -437,14 +446,14 @@ function constitutive_homogenizedC(ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) case (PLASTICITY_DISLOKMC_ID) - constitutive_homogenizedC = plastic_dislokmc_homogenizedC(ipc,ip,el) + constitutive_homogenizedC = plastic_dislokmc_homogenizedC(ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) - constitutive_homogenizedC = plastic_disloucla_homogenizedC(ipc,ip,el) + constitutive_homogenizedC = plastic_disloucla_homogenizedC(ipc,ip,el) case (PLASTICITY_TITANMOD_ID) constitutive_homogenizedC = plastic_titanmod_homogenizedC (ipc,ip,el) case default constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phase (ipc,ip,el)) - + end select end function constitutive_homogenizedC @@ -452,9 +461,9 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) use prec, only: & - pReal + pReal use material, only: & phase_plasticity, & material_phase, & @@ -465,7 +474,8 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) PLASTICITY_dislokmc_ID, & PLASTICITY_disloucla_ID, & PLASTICITY_titanmod_ID, & - PLASTICITY_nonlocal_ID + PLASTICITY_nonlocal_ID, & + PLASTICITY_phenoplus_ID use plastic_titanmod, only: & plastic_titanmod_microstructure use plastic_nonlocal, only: & @@ -476,6 +486,8 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) plastic_dislokmc_microstructure use plastic_disloucla, only: & plastic_disloucla_microstructure + use plastic_phenoplus, only: & + plastic_phenoplus_microstructure implicit none integer(pInt), intent(in) :: & @@ -486,35 +498,39 @@ subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el) Fe, & !< elastic deformation gradient Fp !< plastic deformation gradient integer(pInt) :: & - phase, homog, offset + phase, homog, offset + real(pReal), intent(in), dimension(:,:,:,:) :: & + orientations !< crystal orientation in quaternions phase = material_phase(ipc,ip,el) homog = material_homog( ip,el) offset = thermalMapping(homog)%p(ip,el) select case (phase_plasticity(phase)) - + case (PLASTICITY_DISLOTWIN_ID) call plastic_dislotwin_microstructure(temperature(homog)%p(offset),ipc,ip,el) case (PLASTICITY_DISLOKMC_ID) call plastic_dislokmc_microstructure (temperature(homog)%p(offset),ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) - call plastic_disloucla_microstructure(temperature(homog)%p(offset),ipc,ip,el) + call plastic_disloucla_microstructure(temperature(homog)%p(offset),ipc,ip,el) case (PLASTICITY_TITANMOD_ID) call plastic_titanmod_microstructure (temperature(homog)%p(offset),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) call plastic_nonlocal_microstructure (Fe,Fp,ip,el) + case (PLASTICITY_PHENOPLUS_ID) + call plastic_phenoplus_microstructure(orientations,ipc,ip,el) end select - + end subroutine constitutive_microstructure !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient +!> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v, Fi, ipc, ip, el) use prec, only: & - pReal + pReal use math, only: & math_transpose33, & math_mul33x33, & @@ -530,6 +546,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v PLASTICITY_NONE_ID, & PLASTICITY_J2_ID, & PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOKMC_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -539,6 +556,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v plastic_j2_LpAndItsTangent use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_LpAndItsTangent + use plastic_phenoplus, only: & + plastic_phenoplus_LpAndItsTangent use plastic_dislotwin, only: & plastic_dislotwin_LpAndItsTangent use plastic_dislokmc, only: & @@ -549,7 +568,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v plastic_titanmod_LpAndItsTangent use plastic_nonlocal, only: & plastic_nonlocal_LpAndItsTangent - + implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -564,14 +583,14 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v real(pReal), intent(out), dimension(3,3,3,3) :: & dLp_dTstar3333, & !< derivative of Lp with respect to Tstar (4th-order tensor) dLp_dFi3333 !< derivative of Lp with respect to Fi (4th-order tensor) - real(pReal), dimension(6) :: & + real(pReal), dimension(6) :: & Mstar_v !< Mandel stress work conjugate with Lp real(pReal), dimension(9,9) :: & dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor) real(pReal), dimension(3,3) :: & temp_33 integer(pInt) :: & - i, j, phase, homog, offset + i, j, phase, homog, offset phase = material_phase(ipc,ip,el) homog = material_homog( ip,el) @@ -579,7 +598,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(math_transpose33(Fi),Fi), & math_Mandel6to33(Tstar_v))) select case (phase_plasticity(phase)) - + case (PLASTICITY_NONE_ID) Lp = 0.0_pReal dLp_dMstar = 0.0_pReal @@ -587,6 +606,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v call plastic_j2_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) + case (PLASTICITY_PHENOPLUS_ID) + call plastic_phenoplus_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v,ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) call plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, & temperature(homog)%p(offset), & @@ -602,7 +623,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v case (PLASTICITY_DISLOUCLA_ID) call plastic_disloucla_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, & temperature(homog)%p(offset), & - ipc,ip,el) + ipc,ip,el) case (PLASTICITY_TITANMOD_ID) call plastic_titanmod_LpAndItsTangent(Lp,dLp_dMstar,Mstar_v, & temperature(homog)%p(offset), & @@ -615,21 +636,21 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt dLp_dFi3333(i,j,1:3,1:3) = math_mul33x33(temp_33,math_transpose33(dLp_dTstar3333(i,j,1:3,1:3))) + & math_mul33x33(math_mul33x33(Fi,dLp_dTstar3333(i,j,1:3,1:3)),math_Mandel6to33(Tstar_v)) - enddo; enddo + enddo; enddo temp_33 = math_mul33x33(math_transpose33(Fi),Fi) do i = 1_pInt, 3_pInt; do j = 1_pInt, 3_pInt dLp_dTstar3333(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar3333(i,j,1:3,1:3)) - enddo; enddo - + enddo; enddo + end subroutine constitutive_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient +!> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v, Fi, ipc, ip, el) use prec, only: & - pReal + pReal use math, only: & math_I3, & math_inv33, & @@ -655,7 +676,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v kinematics_vacancy_strain_LiAndItsTangent use kinematics_hydrogen_strain, only: & kinematics_hydrogen_strain_LiAndItsTangent - + implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -673,43 +694,43 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v real(pReal), dimension(3,3) :: & my_Li !< intermediate velocity gradient real(pReal), dimension(3,3,3,3) :: & - my_dLi_dTstar + my_dLi_dTstar real(pReal), dimension(3,3) :: & FiInv, & temp_33 real(pReal) :: & detFi integer(pInt) :: & - i, j, kinematics + i, j, kinematics Li = 0.0_pReal dLi_dTstar3333 = 0.0_pReal dLi_dFi3333 = 0.0_pReal - + do kinematics = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) select case (phase_kinematics(kinematics,material_phase(ipc,ip,el))) case (KINEMATICS_cleavage_opening_ID) call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) - + case (KINEMATICS_slipplane_opening_ID) call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el) - + case (KINEMATICS_thermal_expansion_ID) call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) - + case (KINEMATICS_vacancy_strain_ID) call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) - + case (KINEMATICS_hydrogen_strain_ID) call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el) - + case default my_Li = 0.0_pReal my_dLi_dTstar = 0.0_pReal end select Li = Li + my_Li dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar - enddo + enddo FiInv = math_inv33(Fi) detFi = math_det33(Fi) @@ -719,17 +740,17 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v dLi_dTstar3333(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dTstar3333(1:3,1:3,i,j)),FiInv)*detFi dLi_dFi3333 (1:3,1:3,i,j) = dLi_dFi3333(1:3,1:3,i,j) + Li*FiInv(j,i) dLi_dFi3333 (1:3,i,1:3,j) = dLi_dFi3333(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - enddo; enddo - + enddo; enddo + end subroutine constitutive_LiAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief collects initial intermediate deformation gradient +!> @brief collects initial intermediate deformation gradient !-------------------------------------------------------------------------------------------------- pure function constitutive_initialFi(ipc, ip, el) use prec, only: & - pReal + pReal use math, only: & math_I3, & math_inv33, & @@ -747,7 +768,7 @@ pure function constitutive_initialFi(ipc, ip, el) kinematics_vacancy_strain_initialStrain use kinematics_hydrogen_strain, only: & kinematics_hydrogen_strain_initialStrain - + implicit none integer(pInt), intent(in) :: & ipc, & !< grain number @@ -756,32 +777,32 @@ pure function constitutive_initialFi(ipc, ip, el) real(pReal), dimension(3,3) :: & constitutive_initialFi !< composite initial intermediate deformation gradient integer(pInt) :: & - kinematics + kinematics constitutive_initialFi = math_I3 - + do kinematics = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) !< Warning: small initial strain assumption select case (phase_kinematics(kinematics,material_phase(ipc,ip,el))) case (KINEMATICS_thermal_expansion_ID) constitutive_initialFi = & constitutive_initialFi + kinematics_thermal_expansion_initialStrain(ipc, ip, el) - + case (KINEMATICS_vacancy_strain_ID) constitutive_initialFi = & constitutive_initialFi + kinematics_vacancy_strain_initialStrain(ipc, ip, el) - + case (KINEMATICS_hydrogen_strain_ID) constitutive_initialFi = & constitutive_initialFi + kinematics_hydrogen_strain_initialStrain(ipc, ip, el) - + end select -enddo - +enddo + end function constitutive_initialFi !-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic deformation gradient depending on the selected elastic law (so far no case switch !! because only hooke is implemented !-------------------------------------------------------------------------------------------------- @@ -802,15 +823,15 @@ subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) real(pReal), intent(out), dimension(3,3,3,3) :: & dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - + call constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) - + end subroutine constitutive_TandItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic deformation gradient using hookes law !-------------------------------------------------------------------------------------------------- subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el) @@ -832,7 +853,7 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, damage, & damageMapping, & porosity, & - porosityMapping, & + porosityMapping, & STIFFNESS_DEGRADATION_damage_ID, & STIFFNESS_DEGRADATION_porosity_ID @@ -849,7 +870,7 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, real(pReal), intent(out), dimension(3,3,3,3) :: & dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - + integer(pInt) :: i, j, phase, homog real(pReal), dimension(3,3) :: E real(pReal), dimension(3,3,3,3) :: C @@ -862,30 +883,30 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, case (STIFFNESS_DEGRADATION_damage_ID) C = damage(homog)%p(damageMapping(homog)%p(ip,el))* & damage(homog)%p(damageMapping(homog)%p(ip,el))* & - C - + C + case (STIFFNESS_DEGRADATION_porosity_ID) C = porosity(homog)%p(porosityMapping(homog)%p(ip,el))* & porosity(homog)%p(porosityMapping(homog)%p(ip,el))* & - C + C end select - enddo - + enddo + E = 0.5_pReal*(math_mul33x33(math_transpose33(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration T = math_mul3333xx33(C,math_mul33x33(math_mul33x33(math_transpose33(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration - + dT_dFe = 0.0_pReal forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt) dT_dFe(i,j,1:3,1:3) = & math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),math_transpose33(Fi))) !< dT_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko - dT_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dT_ij/dFi_kl = C_ijln * E_km * Fe_mn + dT_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dT_ij/dFi_kl = C_ijln * E_km * Fe_mn end forall - + end subroutine constitutive_hooke_TandItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfracArray,ipc, ip, el) use prec, only: & @@ -912,6 +933,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra PLASTICITY_none_ID, & PLASTICITY_j2_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & PLASTICITY_disloucla_ID, & @@ -925,6 +947,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra plastic_j2_dotState use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_dotState + use plastic_phenoplus, only: & + plastic_phenoplus_dotState use plastic_dislotwin, only: & plastic_dislotwin_dotState use plastic_dislokmc, only: & @@ -959,15 +983,15 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) integer(pLongInt) :: & - tick, tock, & + tick, tock, & tickrate, & maxticks integer(pInt) :: & phase, homog, offset, mySource - + if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) & call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) - + phase = material_phase(ipc,ip,el) homog = material_homog( ip,el) offset = thermalMapping(homog)%p(ip,el) @@ -976,6 +1000,8 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra call plastic_j2_dotState (Tstar_v,ipc,ip,el) case (PLASTICITY_PHENOPOWERLAW_ID) call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) + case (PLASTICITY_PHENOPLUS_ID) + call plastic_phenoplus_dotState(Tstar_v,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) call plastic_dislotwin_dotState (Tstar_v,temperature(homog)%p(offset), & ipc,ip,el) @@ -992,7 +1018,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,temperature(homog)%p(offset), & subdt,subfracArray,ip,el) end select - + do mySource = 1_pInt, phase_Nsources(phase) select case (phase_source(mySource,phase)) case (SOURCE_damage_anisoBrittle_ID) @@ -1005,7 +1031,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra call source_thermal_externalheat_dotState( ipc, ip, el) end select - enddo + enddo if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) @@ -1040,7 +1066,7 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) PLASTICITY_NONLOCAL_ID, & SOURCE_damage_isoBrittle_ID, & SOURCE_vacancy_irradiation_ID, & - SOURCE_vacancy_thermalfluc_ID + SOURCE_vacancy_thermalfluc_ID use plastic_nonlocal, only: & plastic_nonlocal_deltaState use source_damage_isoBrittle, only: & @@ -1049,20 +1075,20 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) source_vacancy_irradiation_deltaState use source_vacancy_thermalfluc, only: & source_vacancy_thermalfluc_deltaState - + implicit none integer(pInt), intent(in) :: & ipc, & !< grain number ip, & !< integration point number el !< element number real(pReal), intent(in), dimension(6) :: & - Tstar_v !< 2nd Piola-Kirchhoff stress + Tstar_v !< 2nd Piola-Kirchhoff stress real(pReal), intent(in), dimension(3,3) :: & Fe !< elastic deformation gradient integer(pInt) :: & mySource integer(pLongInt) :: & - tick, tock, & + tick, tock, & tickrate, & maxticks @@ -1086,7 +1112,7 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el) call source_vacancy_thermalfluc_deltaState(ipc, ip, el) end select - enddo + enddo if (iand(debug_level(debug_constitutive), debug_levelBasic) /= 0_pInt) then call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) @@ -1106,7 +1132,7 @@ end subroutine constitutive_collectDeltaState !-------------------------------------------------------------------------------------------------- function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) use prec, only: & - pReal + pReal use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -1124,6 +1150,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) PLASTICITY_NONE_ID, & PLASTICITY_J2_ID, & PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_PHENOPLUS_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOKMC_ID, & PLASTICITY_DISLOUCLA_ID, & @@ -1140,6 +1167,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plastic_j2_postResults use plastic_phenopowerlaw, only: & plastic_phenopowerlaw_postResults + use plastic_phenoplus, only: & + plastic_phenoplus_postResults use plastic_dislotwin, only: & plastic_dislotwin_postResults use plastic_dislokmc, only: & @@ -1165,7 +1194,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) ip, & !< integration point number el !< element number real(pReal), dimension(plasticState(material_phase(ipc,ip,el))%sizePostResults + & - sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & + sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & constitutive_postResults real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & FeArray !< elastic deformation gradient @@ -1173,9 +1202,9 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) integer(pInt) :: & startPos, endPos, phase, homog, offset, mySource - + constitutive_postResults = 0.0_pReal - + phase = material_phase(ipc,ip,el) homog = material_homog( ip,el) offset = thermalMapping(homog)%p(ip,el) @@ -1190,6 +1219,9 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) case (PLASTICITY_PHENOPOWERLAW_ID) constitutive_postResults(startPos:endPos) = & plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) + case (PLASTICITY_PHENOPLUS_ID) + constitutive_postResults(startPos:endPos) = & + plastic_phenoplus_postResults(Tstar_v,ipc,ip,el) case (PLASTICITY_DISLOTWIN_ID) constitutive_postResults(startPos:endPos) = & plastic_dislotwin_postResults(Tstar_v,temperature(homog)%p(offset),ipc,ip,el) @@ -1198,7 +1230,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plastic_dislokmc_postResults(Tstar_v,temperature(homog)%p(offset),ipc,ip,el) case (PLASTICITY_DISLOUCLA_ID) constitutive_postResults(startPos:endPos) = & - plastic_disloucla_postResults(Tstar_v,temperature(homog)%p(offset),ipc,ip,el) + plastic_disloucla_postResults(Tstar_v,temperature(homog)%p(offset),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) constitutive_postResults(startPos:endPos) = & plastic_nonlocal_postResults (Tstar_v,FeArray,ip,el) diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 59b847f14..18faf9be7 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -1,9 +1,10 @@ !-------------------------------------------------------------------------------------------------- ! $Id$ !-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH +!> @author Christoph Kords, Max-Planck-Institut für Eisenforschung GmbH +!> @author Chen Zhang, Michigan State University !> @brief crystallite state integration functions and reporting of results !-------------------------------------------------------------------------------------------------- @@ -206,7 +207,7 @@ subroutine crystallite_init tag = '', & line= '' - mainProcess: if (worldrank == 0) then + mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- crystallite init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -406,7 +407,7 @@ subroutine crystallite_init enddo close(FILEUNIT) - endif + endif !-------------------------------------------------------------------------------------------------- ! initialize @@ -419,7 +420,7 @@ subroutine crystallite_init crystallite_F0(1:3,1:3,g,i,e) = math_I3 crystallite_localPlasticity(g,i,e) = phase_localPlasticity(material_phase(g,i,e)) crystallite_Fe(1:3,1:3,g,i,e) = math_inv33(math_mul33x33(crystallite_Fi0(1:3,1:3,g,i,e), & - crystallite_Fp0(1:3,1:3,g,i,e))) ! assuming that euler angles are given in internal strain free configuration + crystallite_Fp0(1:3,1:3,g,i,e))) ! assuming that euler angles are given in internal strain free configuration crystallite_Fp(1:3,1:3,g,i,e) = crystallite_Fp0(1:3,1:3,g,i,e) crystallite_Fi(1:3,1:3,g,i,e) = crystallite_Fi0(1:3,1:3,g,i,e) crystallite_requested(g,i,e) = .true. @@ -437,13 +438,17 @@ subroutine crystallite_init call crystallite_orientations() crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations + !***some debugging statement here + !write(6,*) 'CZ: before crystallite initialization' !$OMP PARALLEL DO PRIVATE(myNgrains) do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do g = 1_pInt,myNgrains + !***dirty way to pass orientation to constitutive module call constitutive_microstructure( & + crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g,i,e) ! update dependent state variables to be consistent with basic states @@ -452,6 +457,9 @@ subroutine crystallite_init enddo !$OMP END PARALLEL DO + !***some debugging statement here + !write(6,*) 'CZ: after crystallite initialization' + call crystallite_stressAndItsTangent(.true.) ! request elastic answers crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback @@ -1132,7 +1140,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_inv33(crystallite_partionedFi0(1:3,1:3,g,i,e))) call constitutive_TandItsTangent(Tstar,dSdFe,dSdFi,Fe_guess,crystallite_partionedFi0(1:3,1:3,g,i,e),g,i,e) crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp), & - math_mul33x33(Tstar,transpose(invFp))) + math_mul33x33(Tstar,transpose(invFp))) endif if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & @@ -1176,7 +1184,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) g,i,e) ! call constitutive law to calculate Li tangent in lattice configuration if (sum(abs(dLidS)) < tol_math_check) then dFidS = 0.0_pReal - else + else temp_33 = math_inv33(crystallite_subFi0(1:3,1:3,g,i,e)) lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal do o=1_pInt,3_pInt; do p=1_pInt,3_pInt @@ -1196,31 +1204,31 @@ subroutine crystallite_stressAndItsTangent(updateJaco) dFidS = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333) endif dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS - endif + endif call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration + crystallite_Fi(1:3,1:3,g,i,e),g,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS temp_33 = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), & - crystallite_invFi(1:3,1:3,g,i,e))) + crystallite_invFi(1:3,1:3,g,i,e))) rhs_3333 = 0.0_pReal forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33) - + temp_3333 = 0.0_pReal temp_33 = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33,dLpdS(1:3,1:3,p,o)), & crystallite_invFi(1:3,1:3,g,i,e)) - + temp_33 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & crystallite_invFp(1:3,1:3,g,i,e)), & math_inv33(crystallite_subFi0(1:3,1:3,g,i,e))) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(temp_33,dLidS(1:3,1:3,p,o)) - + lhs_3333 = crystallite_subdt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + & math_mul3333xx3333(dSdFi,dFidS) @@ -1231,8 +1239,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) dSdF = rhs_3333 else dSdF = math_mul3333xx3333(math_Plain99to3333(temp_99),rhs_3333) - endif - + endif + dFpinvdF = 0.0_pReal temp_3333 = math_mul3333xx3333(dLpdS,dSdF) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & @@ -1240,7 +1248,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)), & math_mul33x33(temp_3333(1:3,1:3,p,o), & crystallite_invFi(1:3,1:3,g,i,e))) - + crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.0_pReal temp_33 = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), & math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), & @@ -1303,14 +1311,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state_backup(:,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state( :,mappingConstitutive(1,g,i,e)) - enddo + enddo plasticState (mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) = & plasticState (mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState_backup(:,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState( :,mappingConstitutive(1,g,i,e)) - enddo + enddo F_backup(1:3,1:3,g,i,e) = crystallite_subF(1:3,1:3,g,i,e) ! ... and kinematics Fp_backup(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) @@ -1348,22 +1356,22 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e); do g = 1,myNgrains plasticState (mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) = & - plasticState (mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) + plasticState (mappingConstitutive(2,g,i,e))%state_backup(:,mappingConstitutive(1,g,i,e)) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state( :,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state_backup(:,mappingConstitutive(1,g,i,e)) - enddo + enddo plasticState (mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = & plasticState (mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState( :,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState_backup(:,mappingConstitutive(1,g,i,e)) - enddo + enddo - crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e) + crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e) crystallite_invFp(1:3,1:3,g,i,e) = InvFp_backup(1:3,1:3,g,i,e) - crystallite_Fi(1:3,1:3,g,i,e) = Fi_backup(1:3,1:3,g,i,e) + crystallite_Fi(1:3,1:3,g,i,e) = Fi_backup(1:3,1:3,g,i,e) crystallite_invFi(1:3,1:3,g,i,e) = InvFi_backup(1:3,1:3,g,i,e) crystallite_Fe(1:3,1:3,g,i,e) = Fe_backup(1:3,1:3,g,i,e) crystallite_Lp(1:3,1:3,g,i,e) = Lp_backup(1:3,1:3,g,i,e) @@ -1383,17 +1391,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state( :,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%subState0(:,mappingConstitutive(1,g,i,e)) - enddo + enddo plasticState (mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = & plasticState (mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState( :,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState_backup(:,mappingConstitutive(1,g,i,e)) - enddo + enddo - crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) - crystallite_Fi(1:3,1:3,g,i,e) = crystallite_subFi0(1:3,1:3,g,i,e) + crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) + crystallite_Fi(1:3,1:3,g,i,e) = crystallite_subFi0(1:3,1:3,g,i,e) crystallite_Fe(1:3,1:3,g,i,e) = crystallite_subFe0(1:3,1:3,g,i,e) crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) crystallite_Li(1:3,1:3,g,i,e) = crystallite_subLi0(1:3,1:3,g,i,e) @@ -1485,14 +1493,14 @@ subroutine crystallite_stressAndItsTangent(updateJaco) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state( :,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%state_backup(:,mappingConstitutive(1,g,i,e)) - enddo + enddo plasticState (mappingConstitutive(2,g,i,e))%dotState( :,mappingConstitutive(1,g,i,e)) = & plasticState (mappingConstitutive(2,g,i,e))%dotState_backup(:,mappingConstitutive(1,g,i,e)) do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState( :,mappingConstitutive(1,g,i,e)) = & sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%dotState_backup(:,mappingConstitutive(1,g,i,e)) - enddo + enddo crystallite_subF(1:3,1:3,g,i,e) = F_backup(1:3,1:3,g,i,e) crystallite_Fp(1:3,1:3,g,i,e) = Fp_backup(1:3,1:3,g,i,e) @@ -1507,7 +1515,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_converged(g,i,e) = convergenceFlag_backup(g,i,e) enddo; enddo enddo elementLooping10 - + deallocate(dPdF_perturbation1) deallocate(dPdF_perturbation2) deallocate(F_backup ) @@ -1596,9 +1604,9 @@ subroutine crystallite_integrateStateRK4() !-------------------------------------------------------------------------------------------------- ! initialize dotState if (.not. singleRun) then - do p = 1_pInt, material_Nphase + do p = 1_pInt, material_Nphase plasticState(p)%RK4dotState = 0.0_pReal - do mySource = 1_pInt, phase_Nsources(p) + do mySource = 1_pInt, phase_Nsources(p) sourceState(p)%p(mySource)%RK4dotState = 0.0_pReal enddo enddo @@ -1606,10 +1614,10 @@ subroutine crystallite_integrateStateRK4() e = eIter(1) i = iIter(1,e) do g = iIter(1,e), iIter(2,e) - plasticState (mappingConstitutive(2,g,i,e))%RK4dotState(:,mappingConstitutive(1,g,i,e)) = 0.0_pReal + plasticState (mappingConstitutive(2,g,i,e))%RK4dotState(:,mappingConstitutive(1,g,i,e)) = 0.0_pReal do mySource = 1_pInt, phase_Nsources(mappingConstitutive(2,g,i,e)) - sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%RK4dotState(:,mappingConstitutive(1,g,i,e)) = 0.0_pReal - enddo + sourceState(mappingConstitutive(2,g,i,e))%p(mySource)%RK4dotState(:,mappingConstitutive(1,g,i,e)) = 0.0_pReal + enddo enddo endif @@ -1630,12 +1638,12 @@ subroutine crystallite_integrateStateRK4() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - c = mappingConstitutive(1,g,i,e) - p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1662,10 +1670,10 @@ subroutine crystallite_integrateStateRK4() p = mappingConstitutive(2,g,i,e) c = mappingConstitutive(1,g,i,e) plasticState(p)%RK4dotState(:,c) = plasticState(p)%RK4dotState(:,c) & - + weight(n)*plasticState(p)%dotState(:,c) + + weight(n)*plasticState(p)%dotState(:,c) do mySource = 1_pInt, phase_Nsources(p) sourceState(p)%p(mySource)%RK4dotState(:,c) = sourceState(p)%p(mySource)%RK4dotState(:,c) & - + weight(n)*sourceState(p)%p(mySource)%dotState(:,c) + + weight(n)*sourceState(p)%p(mySource)%dotState(:,c) enddo endif enddo; enddo; enddo @@ -1729,7 +1737,9 @@ subroutine crystallite_integrateStateRK4() !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) & - call constitutive_microstructure(crystallite_Fe(1:3,1:3,g,i,e), & + !***dirty way to pass orientation information + call constitutive_microstructure(crystallite_orientation, & + crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo @@ -1755,7 +1765,7 @@ subroutine crystallite_integrateStateRK4() ! --- dot state and RK dot state--- - + first3steps: if (n < 4) then !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains @@ -1773,12 +1783,12 @@ subroutine crystallite_integrateStateRK4() !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - p = mappingConstitutive(2,g,i,e) - c = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1946,12 +1956,12 @@ subroutine crystallite_integrateStateRKCK45() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - cc = mappingConstitutive(1,g,i,e) - p = mappingConstitutive(2,g,i,e) + cc = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -1980,7 +1990,7 @@ subroutine crystallite_integrateStateRKCK45() p = mappingConstitutive(2,g,i,e) cc = mappingConstitutive(1,g,i,e) plasticState(p)%RKCK45dotState(stage,:,cc) = plasticState(p)%dotState(:,cc) ! store Runge-Kutta dotState - do mySource = 1_pInt, phase_Nsources(p) + do mySource = 1_pInt, phase_Nsources(p) sourceState(p)%p(mySource)%RKCK45dotState(stage,:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) enddo endif @@ -2054,7 +2064,9 @@ subroutine crystallite_integrateStateRKCK45() !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) & - call constitutive_microstructure(crystallite_Fe(1:3,1:3,g,i,e), & + !***dirty way to pass orientations to constitutive_microstructure + call constitutive_microstructure(crystallite_orientation, & + crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo @@ -2099,12 +2111,12 @@ subroutine crystallite_integrateStateRKCK45() !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - p = mappingConstitutive(2,g,i,e) - cc = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) + cc = mappingConstitutive(1,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,cc))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,cc))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2131,8 +2143,8 @@ subroutine crystallite_integrateStateRKCK45() !$OMP DO PRIVATE(p,cc) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then - p = mappingConstitutive(2,g,i,e) - cc = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) + cc = mappingConstitutive(1,g,i,e) plasticState(p)%RKCK45dotState(6,:,cc) = plasticState (p)%dotState(:,cc) ! store Runge-Kutta dotState do mySource = 1_pInt, phase_Nsources(p) sourceState(p)%p(mySource)%RKCK45dotState(6,:,cc) = sourceState(p)%p(mySource)%dotState(:,cc) ! store Runge-Kutta dotState @@ -2140,7 +2152,7 @@ subroutine crystallite_integrateStateRKCK45() endif enddo; enddo; enddo !$OMP ENDDO - + !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -2157,22 +2169,22 @@ subroutine crystallite_integrateStateRKCK45() sourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e) = & matmul(transpose(sourceState(p)%p(mySource)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),DB) & * crystallite_subdt(g,i,e) - enddo - + enddo + ! --- dot state --- plasticState(p)%dotState(:,cc) = & matmul(transpose(plasticState(p)%RKCK45dotState(1:6,1:mySizePlasticDotState,cc)), B) do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState sourceState(p)%p(mySource)%dotState(:,cc) = & - matmul(transpose(sourceState(p)%p(mySource)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),B) - enddo + matmul(transpose(sourceState(p)%p(mySource)%RKCK45dotState(1:6,1:mySizeSourceDotState,cc)),B) + enddo endif enddo; enddo; enddo !$OMP ENDDO - - ! --- state and update --- - + + ! --- state and update --- + !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -2194,9 +2206,9 @@ subroutine crystallite_integrateStateRKCK45() endif enddo; enddo; enddo !$OMP ENDDO - - ! --- relative residui and state convergence --- - + + ! --- relative residui and state convergence --- + !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,cc,s) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -2206,7 +2218,7 @@ subroutine crystallite_integrateStateRKCK45() forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%state(s,cc)) > 0.0_pReal) & relPlasticStateResiduum(s,g,i,e) = & plasticStateResiduum(s,g,i,e) / plasticState(p)%state(s,cc) - + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState forall (s = 1_pInt:mySizeSourceDotState,abs(sourceState(p)%p(mySource)%state(s,cc)) > 0.0_pReal) & @@ -2220,7 +2232,7 @@ subroutine crystallite_integrateStateRKCK45() rTol_crystalliteState .or. & abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & plasticState(p)%aTolState(1:mySizePlasticDotState)) - do mySource = 1_pInt, phase_Nsources(p) + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState crystallite_todo(g,i,e) = crystallite_todo(g,i,e) .and. & all(abs(relSourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & @@ -2272,7 +2284,9 @@ subroutine crystallite_integrateStateRKCK45() !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) & - call constitutive_microstructure(crystallite_Fe(1:3,1:3,g,i,e), & + !***dirty way to pass orientations to constitutive_microstructure + call constitutive_microstructure(crystallite_orientation, & + crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo @@ -2319,10 +2333,10 @@ subroutine crystallite_integrateStateRKCK45() ! --- nonlocal convergence check --- if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt) & - write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' ! if not requesting Integration of just a single IP + write(6,'(a,i8,a,i2,/)') '<< CRYST >> ', count(crystallite_converged(:,:,:)), ' grains converged' ! if not requesting Integration of just a single IP if ((.not. singleRun) .and. any(.not. crystallite_converged .and. .not. crystallite_localPlasticity)) & ! any non-local not yet converged (or broken)... crystallite_converged = crystallite_converged .and. crystallite_localPlasticity ! ...restart all non-local as not converged - + end subroutine crystallite_integrateStateRKCK45 @@ -2430,12 +2444,12 @@ subroutine crystallite_integrateStateAdaptiveEuler() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - p = mappingConstitutive(2,g,i,e) - c = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2451,7 +2465,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() ! --- STATE UPDATE (EULER INTEGRATION) --- - + !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -2505,7 +2519,9 @@ subroutine crystallite_integrateStateAdaptiveEuler() !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) & - call constitutive_microstructure(crystallite_Fe(1:3,1:3,g,i,e), & + !***dirty way to pass orientations to constitutive_microstructure + call constitutive_microstructure(crystallite_orientation, & + crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo @@ -2555,7 +2571,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2597,8 +2613,8 @@ subroutine crystallite_integrateStateAdaptiveEuler() enddo !$OMP FLUSH(plasticStateResiduum) !$OMP FLUSH(sourceStateResiduum) - - ! --- relative residui --- + + ! --- relative residui --- forall (s = 1_pInt:mySizePlasticDotState, abs(plasticState(p)%dotState(s,c)) > 0.0_pReal) & relPlasticStateResiduum(s,g,i,e) = & plasticStateResiduum(s,g,i,e) / plasticState(p)%dotState(s,c) @@ -2610,8 +2626,8 @@ subroutine crystallite_integrateStateAdaptiveEuler() enddo !$OMP FLUSH(relPlasticStateResiduum) !$OMP FLUSH(relSourceStateResiduum) - -#ifndef _OPENMP + +#ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g)& @@ -2632,7 +2648,7 @@ subroutine crystallite_integrateStateAdaptiveEuler() rTol_crystalliteState .or. & abs(plasticStateResiduum(1:mySizePlasticDotState,g,i,e)) < & plasticState(p)%aTolState(1:mySizePlasticDotState)) - do mySource = 1_pInt, phase_Nsources(p) + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState converged = converged .and. & all(abs(relSourceStateResiduum(1:mySizeSourceDotState,mySource,g,i,e)) < & @@ -2767,12 +2783,12 @@ eIter = FEsolving_execElem(1:2) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - c = mappingConstitutive(1,g,i,e) - p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e) .and. .not. numerics_timeSyncing) then ! if broken non-local... !$OMP CRITICAL (checkTodo) @@ -2788,7 +2804,7 @@ eIter = FEsolving_execElem(1:2) ! --- UPDATE STATE --- - + !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then @@ -2798,7 +2814,7 @@ eIter = FEsolving_execElem(1:2) plasticState(p)%state( 1:mySizePlasticDotState,c) = & plasticState(p)%state( 1:mySizePlasticDotState,c) & + plasticState(p)%dotState(1:mySizePlasticDotState,c) & - * crystallite_subdt(g,i,e) + * crystallite_subdt(g,i,e) do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState sourceState(p)%p(mySource)%state( 1:mySizeSourceDotState,c) = & @@ -2847,7 +2863,9 @@ eIter = FEsolving_execElem(1:2) !$OMP DO do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & - call constitutive_microstructure(crystallite_Fe(1:3,1:3,g,i,e), & + !***dirty way to pass orientations to constitutive_microstructure + call constitutive_microstructure(crystallite_orientation, & + crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g, i, e) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo @@ -2990,7 +3008,7 @@ subroutine crystallite_integrateStateFPI() !-------------------------------------------------------------------------------------------------- ! initialize dotState if (.not. singleRun) then - forall(p = 1_pInt:size(plasticState)) + forall(p = 1_pInt:size(plasticState)) plasticState(p)%previousDotState = 0.0_pReal plasticState(p)%previousDotState2 = 0.0_pReal end forall @@ -3002,13 +3020,13 @@ subroutine crystallite_integrateStateFPI() e = eIter(1) i = iIter(1,e) do g = gIter(1,e), gIter(2,e) - p = mappingConstitutive(2,g,i,e) - c = mappingConstitutive(1,g,i,e) - plasticState(p)%previousDotState (:,c) = 0.0_pReal - plasticState(p)%previousDotState2(:,c) = 0.0_pReal + p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) + plasticState(p)%previousDotState (:,c) = 0.0_pReal + plasticState(p)%previousDotState2(:,c) = 0.0_pReal do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%previousDotState (:,c) = 0.0_pReal - sourceState(p)%p(mySource)%previousDotState2(:,c) = 0.0_pReal + sourceState(p)%p(mySource)%previousDotState (:,c) = 0.0_pReal + sourceState(p)%p(mySource)%previousDotState2(:,c) = 0.0_pReal enddo enddo endif @@ -3032,12 +3050,12 @@ subroutine crystallite_integrateStateFPI() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e)) then - p = mappingConstitutive(2,g,i,e) - c = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) - enddo + enddo if (NaN) then ! NaN occured in any dotState if (.not. crystallite_localPlasticity(g,i,e)) then ! if broken is a non-local... !$OMP CRITICAL (checkTodo) @@ -3052,7 +3070,7 @@ subroutine crystallite_integrateStateFPI() !$OMP ENDDO ! --- UPDATE STATE --- - + !$OMP DO PRIVATE(mySizePlasticDotState,mySizeSourceDotState,p,c) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e)) then @@ -3090,16 +3108,18 @@ subroutine crystallite_integrateStateFPI() !$OMP DO PRIVATE(p,c) do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & - call constitutive_microstructure(crystallite_Fe(1:3,1:3,g,i,e), & + !***dirty way to pass orientations to constitutive_micrsotructure + call constitutive_microstructure(crystallite_orientation, & + crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & g, i, e) ! update dependent state variables to be consistent with basic states - p = mappingConstitutive(2,g,i,e) - c = mappingConstitutive(1,g,i,e) - plasticState(p)%previousDotState2(:,c) = plasticState(p)%previousDotState(:,c) - plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c) + p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) + plasticState(p)%previousDotState2(:,c) = plasticState(p)%previousDotState(:,c) + plasticState(p)%previousDotState (:,c) = plasticState(p)%dotState(:,c) do mySource = 1_pInt, phase_Nsources(p) - sourceState(p)%p(mySource)%previousDotState2(:,c) = sourceState(p)%p(mySource)%previousDotState(:,c) - sourceState(p)%p(mySource)%previousDotState (:,c) = sourceState(p)%p(mySource)%dotState(:,c) + sourceState(p)%p(mySource)%previousDotState2(:,c) = sourceState(p)%p(mySource)%previousDotState(:,c) + sourceState(p)%p(mySource)%previousDotState (:,c) = sourceState(p)%p(mySource)%dotState(:,c) enddo enddo; enddo; enddo !$OMP ENDDO @@ -3145,12 +3165,12 @@ subroutine crystallite_integrateStateFPI() do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains !$OMP FLUSH(crystallite_todo) if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) then - p = mappingConstitutive(2,g,i,e) - c = mappingConstitutive(1,g,i,e) + p = mappingConstitutive(2,g,i,e) + c = mappingConstitutive(1,g,i,e) NaN = any(prec_isNaN(plasticState(p)%dotState(:,c))) do mySource = 1_pInt, phase_Nsources(p) NaN = NaN .or. any(prec_isNaN(sourceState(p)%p(mySource)%dotState(:,c))) - enddo + enddo if (NaN) then ! NaN occured in any dotState crystallite_todo(g,i,e) = .false. ! ... skip me next time if (.not. crystallite_localPlasticity(g,i,e)) then ! if me is non-local... @@ -3207,13 +3227,13 @@ subroutine crystallite_integrateStateFPI() tempPlasticState(1:mySizePlasticDotState) = & plasticState(p)%state(1:mySizePlasticDotState,c) & - plasticStateResiduum(1:mySizePlasticDotState) ! need to copy to local variable, since we cant flush a pointer in openmp - + ! --- store corrected dotState --- (cannot do this before state update, because not sure how to flush pointers in openmp) plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * plasticStateDamper & + plasticState(p)%previousDotState(:,c) & * (1.0_pReal - plasticStateDamper) - + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState dot_prod12 = dot_product( sourceState(p)%p(mySource)%dotState (:,c) & @@ -3252,7 +3272,7 @@ subroutine crystallite_integrateStateFPI() sourceState(p)%p(mySource)%dotState(:,c) * sourceStateDamper & + sourceState(p)%p(mySource)%previousDotState(:,c) & * (1.0_pReal - sourceStateDamper) - enddo + enddo #ifndef _OPENMP if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & @@ -3270,13 +3290,13 @@ subroutine crystallite_integrateStateFPI() plasticState(p)%aTolState(1:mySizePlasticDotState) & .or. abs(plasticStateResiduum(1:mySizePlasticDotState)) < & rTol_crystalliteState * abs(tempPlasticState(1:mySizePlasticDotState))) - do mySource = 1_pInt, phase_Nsources(p) + do mySource = 1_pInt, phase_Nsources(p) mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState converged = converged .and. & all( abs(sourceStateResiduum(1:mySizeSourceDotState,mySource)) < & sourceState(p)%p(mySource)%aTolState(1:mySizeSourceDotState) & .or. abs(sourceStateResiduum(1:mySizeSourceDotState,mySource)) < & - rTol_crystalliteState * abs(tempSourceState(1:mySizeSourceDotState,mySource))) + rTol_crystalliteState * abs(tempSourceState(1:mySizeSourceDotState,mySource))) enddo if (converged) then crystallite_converged(g,i,e) = .true. ! ... converged per definition @@ -3294,7 +3314,7 @@ subroutine crystallite_integrateStateFPI() mySizeSourceDotState = sourceState(p)%p(mySource)%sizeDotState sourceState(p)%p(mySource)%state(1:mySizeSourceDotState,c) = & tempSourceState(1:mySizeSourceDotState,mySource) - enddo + enddo endif enddo; enddo; enddo !$OMP ENDDO @@ -3677,7 +3697,7 @@ logical function crystallite_integrateStress(& residuumLp_old = 0.0_pReal Lpguess_old = Lpguess - LpLoop: do ! inner stress integration loop for consistency with Fi + LpLoop: do ! inner stress integration loop for consistency with Fi NiterationStressLp = NiterationStressLp + 1_pInt loopsExeced: if (NiterationStressLp > nStress) then #ifndef _OPENMP @@ -3687,7 +3707,7 @@ logical function crystallite_integrateStress(& #endif return endif loopsExeced - + !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law B = math_I3 - dt*Lpguess @@ -3889,7 +3909,7 @@ logical function crystallite_integrateStress(& jacoCounterLi = jacoCounterLi + 1_pInt ! increase counter for jaco update Liguess = Liguess + steplengthLi * deltaLi - enddo LiLoop + enddo LiLoop if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (distributionStress) @@ -4034,7 +4054,7 @@ subroutine crystallite_orientations !$OMP PARALLEL DO PRIVATE(myPhase,neighboring_e,neighboring_i,neighboringPhase) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) - myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point) + myPhase = material_phase(1,i,e) ! get my phase (non-local models make no sense with more than one grain per material point) if (plasticState(myPhase)%nonLocal) then ! if nonlocal model ! --- calculate disorientation between me and my neighbor --- @@ -4115,7 +4135,7 @@ function crystallite_postResults(ipc, ip, el) real(pReal), dimension(1+crystallite_sizePostResults(microstructure_crystallite(mesh_element(4,el))) + & 1+plasticState(material_phase(ipc,ip,el))%sizePostResults + & - sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & + sum(sourceState(material_phase(ipc,ip,el))%p(:)%sizePostResults)) :: & crystallite_postResults real(pReal), dimension(3,3) :: & Ee @@ -4129,7 +4149,7 @@ function crystallite_postResults(ipc, ip, el) crystID, & mySize, & n - + crystID = microstructure_crystallite(mesh_element(4,el)) diff --git a/code/debug.f90 b/code/debug.f90 index 19dc8eca9..8aa3fbc5c 100644 --- a/code/debug.f90 +++ b/code/debug.f90 @@ -44,7 +44,7 @@ module debug debug_ABAQUS = 13_pInt integer(pInt), parameter, private :: & debug_MAXNTYPE = debug_ABAQUS !< must be set to the maximum defined debug type - + integer(pInt),protected, dimension(debug_maxNtype+2_pInt), public :: & ! specific ones, and 2 for "all" and "other" debug_level = 0_pInt @@ -60,7 +60,7 @@ module debug debug_cumLpTicks = 0_pLongInt, & !< total cpu ticks spent in LpAndItsTangent debug_cumDeltaStateTicks = 0_pLongInt, & !< total cpu ticks spent in deltaState debug_cumDotStateTicks = 0_pLongInt !< total cpu ticks spent in dotState - + integer(pInt), dimension(2), public :: & debug_stressMaxLocation = 0_pInt, & debug_stressMinLocation = 0_pInt, & @@ -76,19 +76,19 @@ module debug debug_StressLoopLiDistribution, & !< distribution of stress iterations until convergence debug_StressLoopLpDistribution, & !< distribution of stress iterations until convergence debug_StateLoopDistribution !< distribution of state iterations until convergence - + real(pReal), public :: & debug_stressMax = -huge(1.0_pReal), & debug_stressMin = huge(1.0_pReal), & debug_jacobianMax = -huge(1.0_pReal), & debug_jacobianMin = huge(1.0_pReal) - + character(len=64), parameter, private :: & debug_CONFIGFILE = 'debug.config' !< name of configuration file -#ifdef PETSc +#ifdef PETSc character(len=1024), parameter, public :: & - PETSCDEBUG = ' -snes_view -snes_monitor ' + PETSCDEBUG = ' -snes_view -snes_monitor ' #endif public :: debug_init, & debug_reset, & @@ -123,19 +123,19 @@ subroutine debug_init IO_EOF implicit none - integer(pInt), parameter :: FILEUNIT = 300_pInt - + integer(pInt), parameter :: FILEUNIT = 300_pInt + integer(pInt) :: i, what integer(pInt), allocatable, dimension(:) :: chunkPos character(len=65536) :: tag, line - mainProcess: if (worldrank == 0) then + mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- debug init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" endif mainProcess - + if (allocated(debug_StressLoopLpDistribution)) & deallocate(debug_StressLoopLpDistribution) allocate(debug_StressLoopLpDistribution(nStress+1,2)) @@ -160,7 +160,7 @@ subroutine debug_init deallocate(debug_MaterialpointLoopDistribution) allocate(debug_MaterialpointLoopDistribution(nHomog+1)) debug_MaterialpointLoopDistribution = 0_pInt - + !-------------------------------------------------------------------------------------------------- ! try to open the config file @@ -179,7 +179,7 @@ subroutine debug_init case ('grain','g','gr') debug_g = IO_intValue(line,chunkPos,2_pInt) end select - + what = 0_pInt select case(tag) case ('debug') @@ -237,14 +237,14 @@ subroutine debug_init endif enddo close(FILEUNIT) - + do i = 1_pInt, debug_maxNtype if (debug_level(i) == 0) & - debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 2_pInt)) ! fill undefined debug types with levels specified by "other" + debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 2_pInt)) ! fill undefined debug types with levels specified by "other" - debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 1_pInt)) ! fill all debug types with levels specified by "all" + debug_level(i) = ior(debug_level(i), debug_level(debug_MAXNTYPE + 1_pInt)) ! fill all debug types with levels specified by "all" enddo - + if (iand(debug_level(debug_debug),debug_LEVELBASIC) /= 0) & write(6,'(a,/)') ' using values from config file' else fileExists @@ -284,7 +284,7 @@ subroutine debug_init case (debug_ABAQUS) tag = ' ABAQUS FEM solver' end select - + if(debug_level(i) /= 0) then write(6,'(3a)') ' debug level for ', trim(tag), ':' if(iand(debug_level(i),debug_LEVELBASIC) /= 0) write(6,'(a)') ' basic' @@ -305,7 +305,7 @@ subroutine debug_init endif end subroutine debug_init - + !-------------------------------------------------------------------------------------------------- !> @brief resets all debug values @@ -380,7 +380,7 @@ subroutine debug_info write(6,'(a33,1x,f12.6)') 'avg CPU time/microsecs per call :',& real(debug_cumDeltaStateTicks,pReal)*1.0e6_pReal/real(tickrate*debug_cumDeltaStateCalls,pReal) endif - + integral = 0_pInt write(6,'(3/,a)') 'distribution_StressLoopLp : stress stiffness' do j=1_pInt,nStress+1_pInt @@ -394,7 +394,7 @@ subroutine debug_info enddo write(6,'(a15,i10,2(1x,i10))') ' total',integral,sum(debug_StressLoopLpDistribution(:,1)), & sum(debug_StressLoopLpDistribution(:,2)) - + integral = 0_pInt write(6,'(3/,a)') 'distribution_StressLoopLi : stress stiffness' do j=1_pInt,nStress+1_pInt @@ -408,7 +408,7 @@ subroutine debug_info enddo write(6,'(a15,i10,2(1x,i10))') ' total',integral,sum(debug_StressLoopLiDistribution(:,1)), & sum(debug_StressLoopLiDistribution(:,2)) - + integral = 0_pInt write(6,'(2/,a)') 'distribution_CrystalliteStateLoop :' do j=1_pInt,nState+1_pInt @@ -422,7 +422,7 @@ subroutine debug_info enddo write(6,'(a15,i10,2(1x,i10))') ' total',integral,sum(debug_StateLoopDistribution(:,1)), & sum(debug_StateLoopDistribution(:,2)) - + integral = 0_pInt write(6,'(2/,a)') 'distribution_CrystalliteCutbackLoop :' do j=1_pInt,nCryst+1_pInt @@ -435,7 +435,7 @@ subroutine debug_info enddo write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution) endif debugOutputCryst - + debugOutputHomog: if (iand(debug_level(debug_HOMOGENIZATION),debug_LEVELBASIC) /= 0) then integral = 0_pInt write(6,'(2/,a)') 'distribution_MaterialpointStateLoop :' @@ -445,8 +445,8 @@ subroutine debug_info write(6,'(i25,1x,i10)') j,debug_MaterialpointStateLoopDistribution(j) endif enddo - write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution) - + write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointStateLoopDistribution) + integral = 0_pInt write(6,'(2/,a)') 'distribution_MaterialpointCutbackLoop :' do j=1_pInt,nHomog+1_pInt @@ -457,19 +457,19 @@ subroutine debug_info write(6,'(i25,a1,i10)') min(nHomog,j),exceed,debug_MaterialpointLoopDistribution(j) endif enddo - write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution) + write(6,'(a15,i10,1x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution) endif debugOutputHomog - + debugOutputCPFEM: if (iand(debug_level(debug_CPFEM),debug_LEVELBASIC) /= 0) then write(6,'(2/,a,/)') ' Extreme values of returned stress and jacobian' write(6,'(a39)') ' value el ip' write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' stress min :', debug_stressMin, debug_stressMinLocation write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' max :', debug_stressMax, debug_stressMaxLocation write(6,'(a14,1x,e12.3,1x,i6,1x,i4)') ' jacobian min :', debug_jacobianMin, debug_jacobianMinLocation - write(6,'(a14,1x,e12.3,1x,i6,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation + write(6,'(a14,1x,e12.3,1x,i6,1x,i4,/)') ' max :', debug_jacobianMax, debug_jacobianMaxLocation endif debugOutputCPFEM !$OMP END CRITICAL (write2out) - + end subroutine debug_info - + end module debug diff --git a/code/material.f90 b/code/material.f90 index bebfc27e7..a6b428147 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -27,6 +27,7 @@ module material PLASTICITY_none_label = 'none', & PLASTICITY_j2_label = 'j2', & PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', & + PLASTICITY_phenoplus_label = 'phenoplus', & PLASTICITY_dislotwin_label = 'dislotwin', & PLASTICITY_dislokmc_label = 'dislokmc', & PLASTICITY_disloucla_label = 'disloucla', & @@ -63,11 +64,11 @@ module material HYDROGENFLUX_cahnhilliard_label = 'cahnhilliard', & HOMOGENIZATION_none_label = 'none', & HOMOGENIZATION_isostrain_label = 'isostrain', & - HOMOGENIZATION_rgc_label = 'rgc' - - + HOMOGENIZATION_rgc_label = 'rgc' - enum, bind(c) + + + enum, bind(c) enumerator :: ELASTICITY_undefined_ID, & ELASTICITY_hooke_ID end enum @@ -76,6 +77,7 @@ module material PLASTICITY_none_ID, & PLASTICITY_j2_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & PLASTICITY_disloucla_ID, & @@ -146,33 +148,33 @@ module material end enum character(len=*), parameter, public :: & - MATERIAL_configFile = 'material.config', & !< generic name for material configuration file - MATERIAL_localFileExt = 'materialConfig' !< extension of solver job name depending material configuration file - + MATERIAL_configFile = 'material.config', & !< generic name for material configuration file + MATERIAL_localFileExt = 'materialConfig' !< extension of solver job name depending material configuration file + character(len=*), parameter, public :: & MATERIAL_partHomogenization = 'homogenization', & !< keyword for homogenization part MATERIAL_partCrystallite = 'crystallite', & !< keyword for crystallite part MATERIAL_partPhase = 'phase' !< keyword for phase part integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & - phase_elasticity !< elasticity of each phase + phase_elasticity !< elasticity of each phase integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public, protected :: & - phase_plasticity !< plasticity of each phase + phase_plasticity !< plasticity of each phase integer(kind(THERMAL_isothermal_ID)), dimension(:), allocatable, public, protected :: & - thermal_type !< thermal transport model + thermal_type !< thermal transport model integer(kind(DAMAGE_none_ID)), dimension(:), allocatable, public, protected :: & - damage_type !< nonlocal damage model + damage_type !< nonlocal damage model integer(kind(VACANCYFLUX_isoconc_ID)), dimension(:), allocatable, public, protected :: & - vacancyflux_type !< vacancy transport model + vacancyflux_type !< vacancy transport model integer(kind(POROSITY_none_ID)), dimension(:), allocatable, public, protected :: & - porosity_type !< porosity evolution model + porosity_type !< porosity evolution model integer(kind(HYDROGENFLUX_isoconc_ID)), dimension(:), allocatable, public, protected :: & - hydrogenflux_type !< hydrogen transport model + hydrogenflux_type !< hydrogen transport model integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable, public, protected :: & phase_source, & !< active sources mechanisms of each phase - phase_kinematics, & !< active kinematic mechanisms of each phase - phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase + phase_kinematics, & !< active kinematic mechanisms of each phase + phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable, public, protected :: & homogenization_type !< type of each homogenization @@ -188,7 +190,7 @@ module material material_Nhomogenization, & !< number of homogenizations material_Nmicrostructure, & !< number of microstructures material_Ncrystallite !< number of crystallite settings - + integer(pInt), dimension(:), allocatable, public, protected :: & phase_Nsources, & !< number of source mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase @@ -212,16 +214,16 @@ module material microstructure_crystallite !< crystallite setting ID of each microstructure real(pReal), dimension(:), allocatable, public, protected :: & - thermal_initialT, & !< initial temperature per each homogenization + thermal_initialT, & !< initial temperature per each homogenization damage_initialPhi, & !< initial damage per each homogenization - vacancyflux_initialCv, & !< initial vacancy concentration per each homogenization + vacancyflux_initialCv, & !< initial vacancy concentration per each homogenization porosity_initialPhi, & !< initial posority per each homogenization - hydrogenflux_initialCh !< initial hydrogen concentration per each homogenization + hydrogenflux_initialCh !< initial hydrogen concentration per each homogenization integer(pInt), dimension(:,:,:), allocatable, public :: & material_phase !< phase (index) of each grain,IP,element integer(pInt), dimension(:,:), allocatable, public :: & - material_homog !< homogenization (index) of each IP,element + material_homog !< homogenization (index) of each IP,element type(tPlasticState), allocatable, dimension(:), public :: & plasticState type(tSourceState), allocatable, dimension(:), public :: & @@ -236,12 +238,12 @@ module material integer(pInt), dimension(:,:,:), allocatable, public, protected :: & material_texture !< texture (index) of each grain,IP,element - + real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & material_EulerAngles !< initial orientation of each grain,IP,element - + logical, dimension(:), allocatable, public, protected :: & - microstructure_active, & + microstructure_active, & microstructure_elemhomo, & !< flag to indicate homogeneous microstructure distribution over element's IPs phase_localPlasticity !< flags phases with local constitutive law @@ -249,13 +251,13 @@ module material character(len=*), parameter, private :: & MATERIAL_partMicrostructure = 'microstructure', & !< keyword for microstructure part MATERIAL_partTexture = 'texture' !< keyword for texture part - + character(len=64), dimension(:), allocatable, private :: & microstructure_name, & !< name of each microstructure texture_name !< name of each texture - + character(len=256), dimension(:), allocatable, private :: & - texture_ODFfile !< name of each ODF file + texture_ODFfile !< name of each ODF file integer(pInt), private :: & material_Ntexture, & !< number of textures @@ -268,42 +270,42 @@ module material texture_symmetry, & !< number of symmetric orientations per texture texture_Ngauss, & !< number of Gauss components per texture texture_Nfiber !< number of Fiber components per texture - + integer(pInt), dimension(:,:), allocatable, private :: & microstructure_phase, & !< phase IDs of each microstructure microstructure_texture !< texture IDs of each microstructure - + real(pReal), dimension(:,:), allocatable, private :: & microstructure_fraction !< vol fraction of each constituent in microstructure - + real(pReal), dimension(:,:,:), allocatable, private :: & material_volume, & !< volume of each grain,IP,element texture_Gauss, & !< data of each Gauss component texture_Fiber, & !< data of each Fiber component texture_transformation !< transformation for each texture - + logical, dimension(:), allocatable, private :: & homogenization_active - + integer(pInt), dimension(:,:,:,:), allocatable, public, target :: mappingConstitutive integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingCrystallite - integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field - integer(pInt), dimension(:,:), allocatable, public, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field - + integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field + integer(pInt), dimension(:,:), allocatable, public, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field + type(tHomogMapping), allocatable, dimension(:), public :: & thermalMapping, & !< mapping for thermal state/fields damageMapping, & !< mapping for damage state/fields vacancyfluxMapping, & !< mapping for vacancy conc state/fields porosityMapping, & !< mapping for porosity state/fields hydrogenfluxMapping !< mapping for hydrogen conc state/fields - + type(p_vec), allocatable, dimension(:), public :: & temperature, & !< temperature field damage, & !< damage field vacancyConc, & !< vacancy conc field porosity, & !< porosity field hydrogenConc, & !< hydrogen conc field - temperatureRate, & !< temperature change rate field + temperatureRate, & !< temperature change rate field vacancyConcRate, & !< vacancy conc change field hydrogenConcRate !< hydrogen conc change field @@ -313,6 +315,7 @@ module material PLASTICITY_none_ID, & PLASTICITY_J2_ID, & PLASTICITY_phenopowerlaw_ID, & + PLASTICITY_phenoplus_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & PLASTICITY_disloucla_ID, & @@ -367,7 +370,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief parses material configuration file -!> @details figures out if solverJobName.materialConfig is present, if not looks for +!> @details figures out if solverJobName.materialConfig is present, if not looks for !> material.config !-------------------------------------------------------------------------------------------------- subroutine material_init() @@ -387,10 +390,10 @@ subroutine material_init() mesh_NcpElems, & mesh_element, & FE_Nips, & - FE_geomtype + FE_geomtype use numerics, only: & worldrank - + implicit none integer(pInt), parameter :: FILEUNIT = 200_pInt integer(pInt) :: m,c,h, myDebug, myPhase, myHomog @@ -404,8 +407,8 @@ subroutine material_init() integer(pInt), dimension(:), allocatable :: HomogenizationPosition myDebug = debug_level(debug_material) - - mainProcess: if (worldrank == 0) then + + mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- material init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() @@ -430,34 +433,34 @@ subroutine material_init() allocate(sourceState (material_Nphase)) do myPhase = 1,material_Nphase allocate(sourceState(myPhase)%p(phase_Nsources(myPhase))) - enddo - + enddo + allocate(homogState (material_Nhomogenization)) allocate(thermalState (material_Nhomogenization)) allocate(damageState (material_Nhomogenization)) allocate(vacancyfluxState (material_Nhomogenization)) allocate(porosityState (material_Nhomogenization)) allocate(hydrogenfluxState (material_Nhomogenization)) - + allocate(thermalMapping (material_Nhomogenization)) allocate(damageMapping (material_Nhomogenization)) allocate(vacancyfluxMapping (material_Nhomogenization)) allocate(porosityMapping (material_Nhomogenization)) allocate(hydrogenfluxMapping(material_Nhomogenization)) - + allocate(temperature (material_Nhomogenization)) allocate(damage (material_Nhomogenization)) allocate(vacancyConc (material_Nhomogenization)) allocate(porosity (material_Nhomogenization)) allocate(hydrogenConc (material_Nhomogenization)) - + allocate(temperatureRate (material_Nhomogenization)) allocate(vacancyConcRate (material_Nhomogenization)) allocate(hydrogenConcRate (material_Nhomogenization)) - + do m = 1_pInt,material_Nmicrostructure if(microstructure_crystallite(m) < 1_pInt .or. & - microstructure_crystallite(m) > material_Ncrystallite) & + microstructure_crystallite(m) > material_Ncrystallite) & call IO_error(150_pInt,m,ext_msg='crystallite') if(minval(microstructure_phase(1:microstructure_Nconstituents(m),m)) < 1_pInt .or. & maxval(microstructure_phase(1:microstructure_Nconstituents(m),m)) > material_Nphase) & @@ -465,10 +468,10 @@ subroutine material_init() if(minval(microstructure_texture(1:microstructure_Nconstituents(m),m)) < 1_pInt .or. & maxval(microstructure_texture(1:microstructure_Nconstituents(m),m)) > material_Ntexture) & call IO_error(150_pInt,m,ext_msg='texture') - if(microstructure_Nconstituents(m) < 1_pInt) & + if(microstructure_Nconstituents(m) < 1_pInt) & call IO_error(151_pInt,m) enddo - + debugOut: if (iand(myDebug,debug_levelExtensive) /= 0_pInt) then write(6,'(/,a,/)') ' MATERIAL configuration' write(6,'(a32,1x,a16,1x,a6)') 'homogenization ','type ','grains' @@ -491,14 +494,14 @@ subroutine material_init() endif enddo endif debugOut - + call material_populateGrains allocate(mappingConstitutive (2,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt) allocate(mappingHomogenization (2, mesh_maxNips,mesh_NcpElems),source=0_pInt) allocate(mappingCrystallite (2,homogenization_maxNgrains, mesh_NcpElems),source=0_pInt) allocate(mappingHomogenizationConst( mesh_maxNips,mesh_NcpElems),source=1_pInt) - + allocate(ConstitutivePosition (material_Nphase), source=0_pInt) allocate(HomogenizationPosition(material_Nhomogenization),source=0_pInt) allocate(CrystallitePosition (material_Nphase), source=0_pInt) @@ -532,7 +535,7 @@ subroutine material_init() allocate(vacancyConcRate (myHomog)%p(1), source=0.0_pReal) allocate(hydrogenConcRate(myHomog)%p(1), source=0.0_pReal) enddo - + end subroutine material_init @@ -556,25 +559,25 @@ subroutine material_parseHomogenization(fileUnit,myPart) IO_EOF use mesh, only: & mesh_element - + implicit none character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit - - + + integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: Nsections, section, s, p character(len=65536) :: & tag, line logical :: echo - - echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') + + echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') Nsections = IO_countSections(fileUnit,myPart) material_Nhomogenization = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) - + allocate(homogenization_name(Nsections)); homogenization_name = '' - allocate(homogenization_type(Nsections), source=HOMOGENIZATION_undefined_ID) + allocate(homogenization_type(Nsections), source=HOMOGENIZATION_undefined_ID) allocate(thermal_type(Nsections), source=THERMAL_isothermal_ID) allocate(damage_type (Nsections), source=DAMAGE_none_ID) allocate(vacancyflux_type(Nsections), source=VACANCYFLUX_isoconc_ID) @@ -597,7 +600,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) forall (s = 1_pInt:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model? Homogenization view, maximum operations depend on maximum number of homog schemes homogenization_Noutput = IO_countTagInPart(fileUnit,myPart,'(output)',Nsections) - + rewind(fileUnit) line = '' ! to have it initialized section = 0_pInt ! - " - @@ -611,7 +614,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read - exit + exit endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section @@ -628,7 +631,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) homogenization_type(section) = HOMOGENIZATION_NONE_ID homogenization_Ngrains(section) = 1_pInt case(HOMOGENIZATION_ISOSTRAIN_label) - homogenization_type(section) = HOMOGENIZATION_ISOSTRAIN_ID + homogenization_type(section) = HOMOGENIZATION_ISOSTRAIN_ID case(HOMOGENIZATION_RGC_label) homogenization_type(section) = HOMOGENIZATION_RGC_ID case default @@ -639,7 +642,7 @@ subroutine material_parseHomogenization(fileUnit,myPart) case ('thermal') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) case(THERMAL_isothermal_label) - thermal_type(section) = THERMAL_isothermal_ID + thermal_type(section) = THERMAL_isothermal_ID case(THERMAL_adiabatic_label) thermal_type(section) = THERMAL_adiabatic_ID case(THERMAL_conduction_label) @@ -651,19 +654,19 @@ subroutine material_parseHomogenization(fileUnit,myPart) case ('damage') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) case(DAMAGE_NONE_label) - damage_type(section) = DAMAGE_none_ID + damage_type(section) = DAMAGE_none_ID case(DAMAGE_LOCAL_label) - damage_type(section) = DAMAGE_local_ID + damage_type(section) = DAMAGE_local_ID case(DAMAGE_NONLOCAL_label) damage_type(section) = DAMAGE_nonlocal_ID case default call IO_error(500_pInt,ext_msg=trim(IO_stringValue(line,chunkPos,2_pInt))) end select - + case ('vacancyflux') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) case(VACANCYFLUX_isoconc_label) - vacancyflux_type(section) = VACANCYFLUX_isoconc_ID + vacancyflux_type(section) = VACANCYFLUX_isoconc_ID case(VACANCYFLUX_isochempot_label) vacancyflux_type(section) = VACANCYFLUX_isochempot_ID case(VACANCYFLUX_cahnhilliard_label) @@ -675,17 +678,17 @@ subroutine material_parseHomogenization(fileUnit,myPart) case ('porosity') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) case(POROSITY_NONE_label) - porosity_type(section) = POROSITY_none_ID + porosity_type(section) = POROSITY_none_ID case(POROSITY_phasefield_label) - porosity_type(section) = POROSITY_phasefield_ID + porosity_type(section) = POROSITY_phasefield_ID case default call IO_error(500_pInt,ext_msg=trim(IO_stringValue(line,chunkPos,2_pInt))) end select - + case ('hydrogenflux') select case (IO_lc(IO_stringValue(line,chunkPos,2_pInt))) case(HYDROGENFLUX_isoconc_label) - hydrogenflux_type(section) = HYDROGENFLUX_isoconc_ID + hydrogenflux_type(section) = HYDROGENFLUX_isoconc_ID case(HYDROGENFLUX_cahnhilliard_label) hydrogenflux_type(section) = HYDROGENFLUX_cahnhilliard_ID case default @@ -697,19 +700,19 @@ subroutine material_parseHomogenization(fileUnit,myPart) case ('initialtemperature','initialt') thermal_initialT(section) = IO_floatValue(line,chunkPos,2_pInt) - + case ('initialdamage') damage_initialPhi(section) = IO_floatValue(line,chunkPos,2_pInt) - + case ('initialvacancyconc','initialcv') vacancyflux_initialCv(section) = IO_floatValue(line,chunkPos,2_pInt) - + case ('initialporosity') porosity_initialPhi(section) = IO_floatValue(line,chunkPos,2_pInt) - + case ('initialhydrogenconc','initialch') hydrogenflux_initialCh(section) = IO_floatValue(line,chunkPos,2_pInt) - + end select endif enddo @@ -736,12 +739,12 @@ subroutine material_parseMicrostructure(fileUnit,myPart) use mesh, only: & mesh_element, & mesh_NcpElems - + implicit none character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit - - + + integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: Nsections, section, constituent, e, i character(len=65536) :: & @@ -764,7 +767,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) call IO_error(155_pInt,ext_msg='Microstructure in geometry > Sections in material.config') forall (e = 1_pInt:mesh_NcpElems) microstructure_active(mesh_element(4,e)) = .true. ! current microstructure used in model? Elementwise view, maximum N operations for N elements - + microstructure_Nconstituents = IO_countTagInPart(fileUnit,myPart,'(constituent)',Nsections) microstructure_maxNconstituents = maxval(microstructure_Nconstituents) microstructure_elemhomo = IO_spotTagInPart(fileUnit,myPart,'/elementhomogeneous/',Nsections) @@ -772,12 +775,12 @@ subroutine material_parseMicrostructure(fileUnit,myPart) allocate(microstructure_phase (microstructure_maxNconstituents,Nsections),source=0_pInt) allocate(microstructure_texture (microstructure_maxNconstituents,Nsections),source=0_pInt) allocate(microstructure_fraction(microstructure_maxNconstituents,Nsections),source=0.0_pReal) - + rewind(fileUnit) line = '' ! to have it initialized section = 0_pInt ! - " - constituent = 0_pInt ! - " - - + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to line = IO_read(fileUnit) enddo @@ -788,7 +791,7 @@ subroutine material_parseMicrostructure(fileUnit,myPart) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read - exit + exit endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section @@ -840,14 +843,14 @@ subroutine material_parseCrystallite(fileUnit,myPart) implicit none character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit - + integer(pInt) :: Nsections, & section character(len=65536) :: line logical :: echo echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') - + Nsections = IO_countSections(fileUnit,myPart) material_Ncrystallite = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -856,7 +859,7 @@ subroutine material_parseCrystallite(fileUnit,myPart) allocate(crystallite_Noutput(Nsections), source=0_pInt) crystallite_Noutput = IO_countTagInPart(fileUnit,myPart,'(output)',Nsections) - + rewind(fileUnit) line = '' ! to have it initialized section = 0_pInt ! - " - @@ -870,7 +873,7 @@ subroutine material_parseCrystallite(fileUnit,myPart) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read - exit + exit endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section @@ -903,16 +906,16 @@ subroutine material_parsePhase(fileUnit,myPart) implicit none character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit - - + + integer(pInt), allocatable, dimension(:) :: chunkPos - integer(pInt) :: Nsections, section, sourceCtr, kinematicsCtr, stiffDegradationCtr, p + integer(pInt) :: Nsections, section, sourceCtr, kinematicsCtr, stiffDegradationCtr, p character(len=65536) :: & tag,line logical :: echo echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') - + Nsections = IO_countSections(fileUnit,myPart) material_Nphase = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) @@ -933,12 +936,12 @@ subroutine material_parsePhase(fileUnit,myPart) phase_Nkinematics = IO_countTagInPart(fileUnit,myPart,'(kinematics)',Nsections) phase_NstiffnessDegradations = IO_countTagInPart(fileUnit,myPart,'(stiffness_degradation)',Nsections) phase_localPlasticity = .not. IO_spotTagInPart(fileUnit,myPart,'/nonlocal/',Nsections) - + allocate(phase_source(maxval(phase_Nsources),Nsections), source=SOURCE_undefined_ID) allocate(phase_kinematics(maxval(phase_Nkinematics),Nsections), source=KINEMATICS_undefined_ID) allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),Nsections), & source=STIFFNESS_DEGRADATION_undefined_ID) - + rewind(fileUnit) line = '' ! to have it initialized section = 0_pInt ! - " - @@ -952,7 +955,7 @@ subroutine material_parsePhase(fileUnit,myPart) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read - exit + exit endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section @@ -981,6 +984,8 @@ subroutine material_parsePhase(fileUnit,myPart) phase_plasticity(section) = PLASTICITY_J2_ID case (PLASTICITY_PHENOPOWERLAW_label) phase_plasticity(section) = PLASTICITY_PHENOPOWERLAW_ID + case (PLASTICITY_PHENOPLUS_label) + phase_plasticity(section) = PLASTICITY_PHENOPLUS_ID case (PLASTICITY_DISLOTWIN_label) phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID case (PLASTICITY_DISLOKMC_label) @@ -1073,12 +1078,12 @@ subroutine material_parseTexture(fileUnit,myPart) math_sampleRandomOri, & math_I3, & math_inv33 - + implicit none character(len=*), intent(in) :: myPart integer(pInt), intent(in) :: fileUnit - - + + integer(pInt), allocatable, dimension(:) :: chunkPos integer(pInt) :: Nsections, section, gauss, fiber, j character(len=65536) :: tag @@ -1086,14 +1091,14 @@ subroutine material_parseTexture(fileUnit,myPart) logical :: echo echo = IO_globalTagInPart(fileUnit,myPart,'/echo/') - + Nsections = IO_countSections(fileUnit,myPart) material_Ntexture = Nsections if (Nsections < 1_pInt) call IO_error(160_pInt,ext_msg=myPart) allocate(texture_name(Nsections)); texture_name='' allocate(texture_ODFfile(Nsections)); texture_ODFfile='' - allocate(texture_symmetry(Nsections), source=1_pInt) + allocate(texture_symmetry(Nsections), source=1_pInt) allocate(texture_Ngauss(Nsections), source=0_pInt) allocate(texture_Nfiber(Nsections), source=0_pInt) @@ -1104,14 +1109,14 @@ subroutine material_parseTexture(fileUnit,myPart) texture_maxNfiber = maxval(texture_Nfiber) allocate(texture_Gauss (5,texture_maxNgauss,Nsections), source=0.0_pReal) allocate(texture_Fiber (6,texture_maxNfiber,Nsections), source=0.0_pReal) - allocate(texture_transformation(3,3,Nsections), source=0.0_pReal) + allocate(texture_transformation(3,3,Nsections), source=0.0_pReal) texture_transformation = spread(math_I3,3,Nsections) - + rewind(fileUnit) line = '' ! to have in initialized section = 0_pInt ! - " - - gauss = 0_pInt ! - " - - fiber = 0_pInt ! - " - + gauss = 0_pInt ! - " - + fiber = 0_pInt ! - " - do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= myPart) ! wind forward to line = IO_read(fileUnit) enddo @@ -1122,7 +1127,7 @@ subroutine material_parseTexture(fileUnit,myPart) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read - exit + exit endif if (echo) write(6,'(2x,a)') trim(line) ! echo back read lines if (IO_getTag(line,'[',']') /= '') then ! next section @@ -1156,7 +1161,7 @@ subroutine material_parseTexture(fileUnit,myPart) call IO_error(157_pInt,section) end select enddo - + case ('hybridia') textureType texture_ODFfile(section) = IO_stringValue(line,chunkPos,2_pInt) @@ -1170,7 +1175,7 @@ subroutine material_parseTexture(fileUnit,myPart) case default texture_symmetry(section) = 1_pInt end select - + case ('(random)') textureType gauss = gauss + 1_pInt texture_Gauss(1:3,gauss,section) = math_sampleRandomOri() @@ -1258,7 +1263,7 @@ subroutine material_populateGrains debug_level, & debug_material, & debug_levelBasic - + implicit none integer(pInt), dimension (:,:), allocatable :: Ngrains integer(pInt), dimension (microstructure_maxNconstituents) :: & @@ -1280,22 +1285,22 @@ subroutine material_populateGrains type(p_intvec), dimension (:,:), allocatable :: elemsOfHomogMicro ! lists element number in homog, micro array myDebug = debug_level(debug_material) - + allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0.0_pReal) allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(material_homog(mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(material_texture(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt) allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0.0_pReal) - + allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure), source=0_pInt) allocate(Nelems(material_Nhomogenization,material_Nmicrostructure), source=0_pInt) - + ! populating homogenization schemes in each !-------------------------------------------------------------------------------------------------- do e = 1_pInt, mesh_NcpElems material_homog(1_pInt:FE_Nips(FE_geomtype(mesh_element(2,e))),e) = mesh_element(3,e) enddo - + !-------------------------------------------------------------------------------------------------- ! precounting of elements for each homog/micro pair do e = 1_pInt, mesh_NcpElems @@ -1333,12 +1338,12 @@ subroutine material_populateGrains Nelems(homog,micro) = Nelems(homog,micro) + 1_pInt ! total element count elemsOfHomogMicro(homog,micro)%p(Nelems(homog,micro)) = e ! remember elements active in this homog/micro pair enddo elementLooping - + allocate(volumeOfGrain(maxval(Ngrains)), source=0.0_pReal) ! reserve memory for maximum case allocate(phaseOfGrain(maxval(Ngrains)), source=0_pInt) ! reserve memory for maximum case allocate(textureOfGrain(maxval(Ngrains)), source=0_pInt) ! reserve memory for maximum case allocate(orientationOfGrain(3,maxval(Ngrains)),source=0.0_pReal) ! reserve memory for maximum case - + if (iand(myDebug,debug_levelBasic) /= 0_pInt) then !$OMP CRITICAL (write2out) write(6,'(/,a/)') ' MATERIAL grain population' @@ -1512,11 +1517,11 @@ subroutine material_populateGrains orientationOfGrain(1:3,grain+t) = orientationOfGrain(1:3,grain+j) orientationOfGrain(1:3,grain+j) = orientation enddo - + enddo texture !< @todo calc fraction after weighing with volumePerGrain, exchange in MC steps to improve result (humbug at the moment) - + !-------------------------------------------------------------------------------------------------- ! distribute grains of all constituents as accurately as possible to given constituent fractions @@ -1579,7 +1584,7 @@ subroutine material_populateGrains endif ! active homog,micro pair enddo enddo - + deallocate(volumeOfGrain) deallocate(phaseOfGrain) deallocate(textureOfGrain) @@ -1600,7 +1605,7 @@ integer(pInt) pure function material_NconstituentsPhase(matID) implicit none integer(pInt), intent(in) :: matID - + material_NconstituentsPhase = count(microstructure_phase == matID) end function #endif diff --git a/code/plastic_phenoplus.f90 b/code/plastic_phenoplus.f90 new file mode 100644 index 000000000..119248e9d --- /dev/null +++ b/code/plastic_phenoplus.f90 @@ -0,0 +1,1356 @@ +!-------------------------------------------------------------------------------------------------- +! $Id: plastic_phenopowerlaw.f90 4457 2015-09-08 19:44:32Z MPIE\m.diehl $ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Chen Zhang, Michigan State University +!> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw +!! fitting +!-------------------------------------------------------------------------------------------------- +module plastic_phenoplus + use prec, only: & + pReal,& + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_phenoplus_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + plastic_phenoplus_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + plastic_phenoplus_output !< name of each post result output + + integer(pInt), dimension(:), allocatable, target, public :: & + plastic_phenoplus_Noutput !< number of outputs per instance of this constitution + + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_phenoplus_totalNslip, & !< no. of slip system used in simulation + plastic_phenoplus_totalNtwin, & !< no. of twin system used in simulation + plastic_phenoplus_totalNtrans !< no. of trans system used in simulation + + integer(pInt), dimension(:,:), allocatable, private :: & + plastic_phenoplus_Nslip, & !< active number of slip systems per family (input parameter, per family) + plastic_phenoplus_Ntwin, & !< active number of twin systems per family (input parameter, per family) + plastic_phenoplus_Ntrans !< active number of trans systems per family (input parameter, per family) + + real(pReal), dimension(:), allocatable, private :: & + plastic_phenoplus_gdot0_slip, & !< reference shear strain rate for slip (input parameter) + plastic_phenoplus_gdot0_twin, & !< reference shear strain rate for twin (input parameter) + plastic_phenoplus_n_slip, & !< stress exponent for slip (input parameter) + plastic_phenoplus_n_twin, & !< stress exponent for twin (input parameter) + plastic_phenoplus_spr, & !< push-up factor for slip saturation due to twinning + plastic_phenoplus_twinB, & + plastic_phenoplus_twinC, & + plastic_phenoplus_twinD, & + plastic_phenoplus_twinE, & + plastic_phenoplus_h0_SlipSlip, & !< reference hardening slip - slip (input parameter) + plastic_phenoplus_h0_TwinSlip, & !< reference hardening twin - slip (input parameter) + plastic_phenoplus_h0_TwinTwin, & !< reference hardening twin - twin (input parameter) + plastic_phenoplus_a_slip, & + plastic_phenoplus_aTolResistance, & + plastic_phenoplus_aTolShear, & + plastic_phenoplus_aTolTwinfrac, & + plastic_phenoplus_aTolTransfrac, & + plastic_phenoplus_Cnuc, & !< coefficient for strain-induced martensite nucleation + plastic_phenoplus_Cdwp, & !< coefficient for double well potential + plastic_phenoplus_Cgro, & !< coefficient for stress-assisted martensite growth + plastic_phenoplus_deltaG, & !< free energy difference between austensite and martensite [MPa] + plastic_phenoplus_kappa_max !< capped kappa for each slip system + + real(pReal), dimension(:,:), allocatable, private :: & + plastic_phenoplus_tau0_slip, & !< initial critical shear stress for slip (input parameter, per family) + plastic_phenoplus_tau0_twin, & !< initial critical shear stress for twin (input parameter, per family) + plastic_phenoplus_tausat_slip, & !< maximum critical shear stress for slip (input parameter, per family) + plastic_phenoplus_nonSchmidCoeff, & + + plastic_phenoplus_interaction_SlipSlip, & !< interaction factors slip - slip (input parameter) + plastic_phenoplus_interaction_SlipTwin, & !< interaction factors slip - twin (input parameter) + plastic_phenoplus_interaction_TwinSlip, & !< interaction factors twin - slip (input parameter) + plastic_phenoplus_interaction_TwinTwin !< interaction factors twin - twin (input parameter) + + real(pReal), dimension(:,:,:), allocatable, private :: & + plastic_phenoplus_hardeningMatrix_SlipSlip, & + plastic_phenoplus_hardeningMatrix_SlipTwin, & + plastic_phenoplus_hardeningMatrix_TwinSlip, & + plastic_phenoplus_hardeningMatrix_TwinTwin + + enum, bind(c) + enumerator :: undefined_ID, & + resistance_slip_ID, & + accumulatedshear_slip_ID, & + shearrate_slip_ID, & + resolvedstress_slip_ID, & + kappa_slip_ID, & + totalshear_ID, & + resistance_twin_ID, & + accumulatedshear_twin_ID, & + shearrate_twin_ID, & + resolvedstress_twin_ID, & + totalvolfrac_twin_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + plastic_phenoplus_outputID !< ID of each post result output + + public :: & + plastic_phenoplus_init, & + plastic_phenoplus_microstructure, & + plastic_phenoplus_LpAndItsTangent, & + plastic_phenoplus_dotState, & + plastic_phenoplus_postResults + private :: & + plastic_phenoplus_aTolState, & + plastic_phenoplus_stateInit + + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine plastic_phenoplus_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level, & + debug_constitutive,& + debug_levelBasic + use math, only: & + math_Mandel3333to66, & + math_Voigt66to3333 + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + phase_plasticity, & + phase_plasticityInstance, & + phase_Noutput, & + PLASTICITY_PHENOPLUS_label, & + PLASTICITY_PHENOPLUS_ID, & + material_phase, & + plasticState, & + MATERIAL_partPhase + use lattice + use numerics,only: & + analyticJaco, & + worldrank, & + numerics_integrator + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), allocatable, dimension(:) :: chunkPos + integer(pInt) :: & + maxNinstance, & + instance,phase,j,k, f,o, & + Nchunks_SlipSlip = 0_pInt, Nchunks_SlipTwin = 0_pInt, & + Nchunks_TwinSlip = 0_pInt, Nchunks_TwinTwin = 0_pInt, & + Nchunks_SlipFamilies = 0_pInt, Nchunks_TwinFamilies = 0_pInt, & + Nchunks_TransFamilies = 0_pInt, Nchunks_nonSchmid = 0_pInt, & + NipcMyPhase, & + offset_slip, index_myFamily, index_otherFamily, & + mySize=0_pInt,sizeState,sizeDotState, sizeDeltaState + character(len=65536) :: & + tag = '', & + line = '' + real(pReal), dimension(:), allocatable :: tempPerSlip + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPLUS_label//' init -+>>>' + write(6,'(a)') ' $Id: plastic_phenoplus.f90 4457 2015-09-08 19:44:32Z MPIE\m.diehl $' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPLUS_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(plastic_phenoplus_sizePostResults(maxNinstance), source=0_pInt) + allocate(plastic_phenoplus_sizePostResult(maxval(phase_Noutput),maxNinstance), & + source=0_pInt) + allocate(plastic_phenoplus_output(maxval(phase_Noutput),maxNinstance)) + plastic_phenoplus_output = '' + allocate(plastic_phenoplus_outputID(maxval(phase_Noutput),maxNinstance),source=undefined_ID) + allocate(plastic_phenoplus_Noutput(maxNinstance), source=0_pInt) + allocate(plastic_phenoplus_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) + allocate(plastic_phenoplus_Ntwin(lattice_maxNtwinFamily,maxNinstance), source=0_pInt) + allocate(plastic_phenoplus_Ntrans(lattice_maxNtransFamily,maxNinstance),source=0_pInt) + allocate(plastic_phenoplus_totalNslip(maxNinstance), source=0_pInt) + allocate(plastic_phenoplus_totalNtwin(maxNinstance), source=0_pInt) + allocate(plastic_phenoplus_totalNtrans(maxNinstance), source=0_pInt) + allocate(plastic_phenoplus_gdot0_slip(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_n_slip(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_tau0_slip(lattice_maxNslipFamily,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_tausat_slip(lattice_maxNslipFamily,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_gdot0_twin(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_n_twin(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_tau0_twin(lattice_maxNtwinFamily,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_spr(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_twinB(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_twinC(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_twinD(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_twinE(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_h0_SlipSlip(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_h0_TwinSlip(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_h0_TwinTwin(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_interaction_SlipSlip(lattice_maxNinteraction,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_interaction_SlipTwin(lattice_maxNinteraction,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_interaction_TwinSlip(lattice_maxNinteraction,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_interaction_TwinTwin(lattice_maxNinteraction,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_a_slip(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_aTolResistance(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_aTolShear(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_aTolTwinfrac(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_aTolTransfrac(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstance), & + source=0.0_pReal) + allocate(plastic_phenoplus_Cnuc(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_Cdwp(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_Cgro(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_deltaG(maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_kappa_max(maxNinstance), source=0.0_pReal) + + rewind(fileUnit) + phase = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next phase + phase = phase + 1_pInt ! advance phase section counter + if (phase_plasticity(phase) == PLASTICITY_PHENOPLUS_ID) then + Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) ! maximum number of slip families according to lattice type of current phase + Nchunks_TwinFamilies = count(lattice_NtwinSystem(:,phase) > 0_pInt) ! maximum number of twin families according to lattice type of current phase + Nchunks_TransFamilies = count(lattice_NtransSystem(:,phase) > 0_pInt) ! maximum number of trans families according to lattice type of current phase + Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) + Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase)) + Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase)) + Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase)) + Nchunks_nonSchmid = lattice_NnonSchmid(phase) + if(allocated(tempPerSlip)) deallocate(tempPerSlip) + allocate(tempPerSlip(Nchunks_SlipFamilies)) + endif + cycle ! skip to next line + endif + if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_PHENOPLUS_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran + instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase + chunkPos = IO_stringPos(line) + tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt))) + case ('resistance_slip') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = resistance_slip_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('accumulatedshear_slip','accumulated_shear_slip') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = accumulatedshear_slip_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('shearrate_slip') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = shearrate_slip_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('resolvedstress_slip') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = resolvedstress_slip_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('kappa_slip') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = kappa_slip_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('totalshear') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = totalshear_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('resistance_twin') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = resistance_twin_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('accumulatedshear_twin','accumulated_shear_twin') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = accumulatedshear_twin_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('shearrate_twin') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = shearrate_twin_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('resolvedstress_twin') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = resolvedstress_twin_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case ('totalvolfrac_twin') + plastic_phenoplus_Noutput(instance) = plastic_phenoplus_Noutput(instance) + 1_pInt + plastic_phenoplus_outputID(plastic_phenoplus_Noutput(instance),instance) = totalvolfrac_twin_ID + plastic_phenoplus_output(plastic_phenoplus_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,chunkPos,2_pInt)) + case default + + end select +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of slip families + case ('nslip') + if (chunkPos(1) < Nchunks_SlipFamilies + 1_pInt) & + call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + if (chunkPos(1) > Nchunks_SlipFamilies + 1_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + Nchunks_SlipFamilies = chunkPos(1) - 1_pInt ! user specified number of (possibly) active slip families (e.g. 6 0 6 --> 3) + do j = 1_pInt, Nchunks_SlipFamilies + plastic_phenoplus_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) + enddo + case ('tausat_slip','tau0_slip') + tempPerSlip = 0.0_pReal + do j = 1_pInt, Nchunks_SlipFamilies + if (plastic_phenoplus_Nslip(j,instance) > 0_pInt) & + tempPerSlip(j) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + select case(tag) + case ('tausat_slip') + plastic_phenoplus_tausat_slip(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('tau0_slip') + plastic_phenoplus_tau0_slip(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + end select +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of twin families + case ('ntwin') + if (chunkPos(1) < Nchunks_TwinFamilies + 1_pInt) & + call IO_warning(51_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + if (chunkPos(1) > Nchunks_TwinFamilies + 1_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + Nchunks_TwinFamilies = chunkPos(1) - 1_pInt + do j = 1_pInt, Nchunks_TwinFamilies + plastic_phenoplus_Ntwin(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) + enddo + case ('tau0_twin') + do j = 1_pInt, Nchunks_TwinFamilies + if (plastic_phenoplus_Ntwin(j,instance) > 0_pInt) & + plastic_phenoplus_tau0_twin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of transformation families + case ('ntrans') + if (chunkPos(1) < Nchunks_TransFamilies + 1_pInt) & + call IO_warning(51_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + if (chunkPos(1) > Nchunks_TransFamilies + 1_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + Nchunks_TransFamilies = chunkPos(1) - 1_pInt + do j = 1_pInt, Nchunks_TransFamilies + plastic_phenoplus_Ntrans(j,instance) = IO_intValue(line,chunkPos,1_pInt+j) + enddo +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of interactions + case ('interaction_slipslip') + if (chunkPos(1) < 1_pInt + Nchunks_SlipSlip) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + do j = 1_pInt, Nchunks_SlipSlip + plastic_phenoplus_interaction_SlipSlip(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + case ('interaction_sliptwin') + if (chunkPos(1) < 1_pInt + Nchunks_SlipTwin) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + do j = 1_pInt, Nchunks_SlipTwin + plastic_phenoplus_interaction_SlipTwin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + case ('interaction_twinslip') + if (chunkPos(1) < 1_pInt + Nchunks_TwinSlip) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + do j = 1_pInt, Nchunks_TwinSlip + plastic_phenoplus_interaction_TwinSlip(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + case ('interaction_twintwin') + if (chunkPos(1) < 1_pInt + Nchunks_TwinTwin) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + do j = 1_pInt, Nchunks_TwinTwin + plastic_phenoplus_interaction_TwinTwin(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo + case ('nonschmid_coefficients') + if (chunkPos(1) < 1_pInt + Nchunks_nonSchmid) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPLUS_label//')') + do j = 1_pInt,Nchunks_nonSchmid + plastic_phenoplus_nonSchmidCoeff(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j) + enddo +!-------------------------------------------------------------------------------------------------- +! parameters independent of number of slip/twin systems + case ('gdot0_slip') + plastic_phenoplus_gdot0_slip(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('n_slip') + plastic_phenoplus_n_slip(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('a_slip', 'w0_slip') + plastic_phenoplus_a_slip(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('gdot0_twin') + plastic_phenoplus_gdot0_twin(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('n_twin') + plastic_phenoplus_n_twin(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('s_pr') + plastic_phenoplus_spr(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('twin_b') + plastic_phenoplus_twinB(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('twin_c') + plastic_phenoplus_twinC(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('twin_d') + plastic_phenoplus_twinD(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('twin_e') + plastic_phenoplus_twinE(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('h0_slipslip') + plastic_phenoplus_h0_SlipSlip(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('h0_twinslip') + plastic_phenoplus_h0_TwinSlip(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('h0_twintwin') + plastic_phenoplus_h0_TwinTwin(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('atol_resistance') + plastic_phenoplus_aTolResistance(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('atol_shear') + plastic_phenoplus_aTolShear(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('atol_twinfrac') + plastic_phenoplus_aTolTwinfrac(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('atol_transfrac') + plastic_phenoplus_aTolTransfrac(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('kappa_max') + plastic_phenoplus_kappa_max(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('cnuc') + plastic_phenoplus_Cnuc(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('cdwp') + plastic_phenoplus_Cdwp(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('cgro') + plastic_phenoplus_Cgro(instance) = IO_floatValue(line,chunkPos,2_pInt) + case ('deltag') + plastic_phenoplus_deltaG(instance) = IO_floatValue(line,chunkPos,2_pInt) + case default + + end select + endif; endif + enddo parsingFile + + sanityChecks: do phase = 1_pInt, size(phase_plasticity) + myPhase: if (phase_plasticity(phase) == PLASTICITY_phenoplus_ID) then + instance = phase_plasticityInstance(phase) + plastic_phenoplus_Nslip(1:lattice_maxNslipFamily,instance) = & + min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested + plastic_phenoplus_Nslip(1:lattice_maxNslipFamily,instance)) + plastic_phenoplus_Ntwin(1:lattice_maxNtwinFamily,instance) = & + min(lattice_NtwinSystem(1:lattice_maxNtwinFamily,phase),& ! limit active twin systems per family to min of available and requested + plastic_phenoplus_Ntwin(:,instance)) + plastic_phenoplus_totalNslip(instance) = sum(plastic_phenoplus_Nslip(:,instance)) ! how many slip systems altogether + plastic_phenoplus_totalNtwin(instance) = sum(plastic_phenoplus_Ntwin(:,instance)) ! how many twin systems altogether + plastic_phenoplus_totalNtrans(instance) = sum(plastic_phenoplus_Ntrans(:,instance)) ! how many trans systems altogether + + if (any(plastic_phenoplus_tau0_slip(:,instance) < 0.0_pReal .and. & + plastic_phenoplus_Nslip(:,instance) > 0)) & + call IO_error(211_pInt,el=instance,ext_msg='tau0_slip ('//PLASTICITY_PHENOPLUS_label//')') + if (plastic_phenoplus_gdot0_slip(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='gdot0_slip ('//PLASTICITY_PHENOPLUS_label//')') + if (plastic_phenoplus_n_slip(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='n_slip ('//PLASTICITY_PHENOPLUS_label//')') + if (any(plastic_phenoplus_tausat_slip(:,instance) <= 0.0_pReal .and. & + plastic_phenoplus_Nslip(:,instance) > 0)) & + call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPLUS_label//')') + if (any(abs(plastic_phenoplus_a_slip(instance)) <= tiny(0.0_pReal) .and. & + plastic_phenoplus_Nslip(:,instance) > 0)) & + call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPLUS_label//')') + if (any(plastic_phenoplus_tau0_twin(:,instance) < 0.0_pReal .and. & + plastic_phenoplus_Ntwin(:,instance) > 0)) & + call IO_error(211_pInt,el=instance,ext_msg='tau0_twin ('//PLASTICITY_PHENOPLUS_label//')') + if ( plastic_phenoplus_gdot0_twin(instance) <= 0.0_pReal .and. & + any(plastic_phenoplus_Ntwin(:,instance) > 0)) & + call IO_error(211_pInt,el=instance,ext_msg='gdot0_twin ('//PLASTICITY_PHENOPLUS_label//')') + if ( plastic_phenoplus_n_twin(instance) <= 0.0_pReal .and. & + any(plastic_phenoplus_Ntwin(:,instance) > 0)) & + call IO_error(211_pInt,el=instance,ext_msg='n_twin ('//PLASTICITY_PHENOPLUS_label//')') + if (plastic_phenoplus_aTolResistance(instance) <= 0.0_pReal) & + plastic_phenoplus_aTolResistance(instance) = 1.0_pReal ! default absolute tolerance 1 Pa + if (plastic_phenoplus_aTolShear(instance) <= 0.0_pReal) & + plastic_phenoplus_aTolShear(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6 + if (plastic_phenoplus_aTolTwinfrac(instance) <= 0.0_pReal) & + plastic_phenoplus_aTolTwinfrac(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6 + if (plastic_phenoplus_aTolTransfrac(instance) <= 0.0_pReal) & + plastic_phenoplus_aTolTransfrac(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6 + endif myPhase + enddo sanityChecks + +!-------------------------------------------------------------------------------------------------- +! allocation of variables whose size depends on the total number of active slip systems + allocate(plastic_phenoplus_hardeningMatrix_SlipSlip(maxval(plastic_phenoplus_totalNslip),& ! slip resistance from slip activity + maxval(plastic_phenoplus_totalNslip),& + maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_hardeningMatrix_SlipTwin(maxval(plastic_phenoplus_totalNslip),& ! slip resistance from twin activity + maxval(plastic_phenoplus_totalNtwin),& + maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_hardeningMatrix_TwinSlip(maxval(plastic_phenoplus_totalNtwin),& ! twin resistance from slip activity + maxval(plastic_phenoplus_totalNslip),& + maxNinstance), source=0.0_pReal) + allocate(plastic_phenoplus_hardeningMatrix_TwinTwin(maxval(plastic_phenoplus_totalNtwin),& ! twin resistance from twin activity + maxval(plastic_phenoplus_totalNtwin),& + maxNinstance), source=0.0_pReal) + + initializeInstances: do phase = 1_pInt, size(phase_plasticity) ! loop through all phases in material.config + myPhase2: if (phase_plasticity(phase) == PLASTICITY_phenoplus_ID) then ! only consider my phase + NipcMyPhase = count(material_phase == phase) ! number of IPCs containing my phase + instance = phase_plasticityInstance(phase) ! which instance of my phase + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputsLoop: do o = 1_pInt,plastic_phenoplus_Noutput(instance) + select case(plastic_phenoplus_outputID(o,instance)) + case(resistance_slip_ID, & + shearrate_slip_ID, & + accumulatedshear_slip_ID, & + resolvedstress_slip_ID, & + kappa_slip_ID & + ) + mySize = plastic_phenoplus_totalNslip(instance) + case(resistance_twin_ID, & + shearrate_twin_ID, & + accumulatedshear_twin_ID, & + resolvedstress_twin_ID & + ) + mySize = plastic_phenoplus_totalNtwin(instance) + case(totalshear_ID, & + totalvolfrac_twin_ID & + ) + mySize = 1_pInt + case default + end select + + outputFound: if (mySize > 0_pInt) then + plastic_phenoplus_sizePostResult(o,instance) = mySize + plastic_phenoplus_sizePostResults(instance) = plastic_phenoplus_sizePostResults(instance) + mySize + endif outputFound + enddo outputsLoop +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + sizeState = plastic_phenoplus_totalNslip(instance) & ! s_slip + + plastic_phenoplus_totalNtwin(instance) & ! s_twin + + 2_pInt & ! sum(gamma) + sum(f) + + plastic_phenoplus_totalNslip(instance) & ! accshear_slip + + plastic_phenoplus_totalNtwin(instance) & ! accshear_twin + + plastic_phenoplus_totalNslip(instance) ! kappa + + !sizeDotState = sizeState ! same as sizeState + !QUICK FIX: the dotState cannot have redundancy, which could cause unknown error + ! explicitly specify the size of the dotState to avoid this potential + ! memory leak issue. + sizeDotState = plastic_phenoplus_totalNslip(instance) & ! s_slip + + plastic_phenoplus_totalNtwin(instance) & ! s_twin + + 2_pInt & ! sum(gamma) + sum(f) + + plastic_phenoplus_totalNslip(instance) & ! accshear_slip + + plastic_phenoplus_totalNtwin(instance) ! accshear_twin + + sizeDeltaState = 0_pInt + plasticState(phase)%sizeState = sizeState + plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%sizeDeltaState = sizeDeltaState + plasticState(phase)%sizePostResults = plastic_phenoplus_sizePostResults(instance) + plasticState(phase)%nSlip =plastic_phenoplus_totalNslip(instance) + plasticState(phase)%nTwin =plastic_phenoplus_totalNtwin(instance) + plasticState(phase)%nTrans=plastic_phenoplus_totalNtrans(instance) + allocate(plasticState(phase)%aTolState ( sizeState), source=0.0_pReal) + allocate(plasticState(phase)%state0 ( sizeState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%partionedState0 ( sizeState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%subState0 ( sizeState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%state ( sizeState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDeltaState,NipcMyPhase), source=0.0_pReal) + if (.not. analyticJaco) then + allocate(plasticState(phase)%state_backup ( sizeState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%dotState_backup (sizeDotState,NipcMyPhase),source=0.0_pReal) + endif + if (any(numerics_integrator == 1_pInt)) then + allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal) + allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) + + offset_slip = plasticState(phase)%nSlip+plasticState(phase)%nTwin+2_pInt + plasticState(phase)%slipRate => & + plasticState(phase)%dotState(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase) + plasticState(phase)%accumulatedSlip => & + plasticState(phase)%state(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase) + + do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X + index_myFamily = sum(plastic_phenoplus_Nslip(1:f-1_pInt,instance)) + do j = 1_pInt,plastic_phenoplus_Nslip(f,instance) ! loop over (active) systems in my family (slip) + do o = 1_pInt,lattice_maxNslipFamily + index_otherFamily = sum(plastic_phenoplus_Nslip(1:o-1_pInt,instance)) + do k = 1_pInt,plastic_phenoplus_Nslip(o,instance) ! loop over (active) systems in other family (slip) + plastic_phenoplus_hardeningMatrix_SlipSlip(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_phenoplus_interaction_SlipSlip(lattice_interactionSlipSlip( & + sum(lattice_NslipSystem(1:f-1,phase))+j, & + sum(lattice_NslipSystem(1:o-1,phase))+k, & + phase), instance ) + enddo; enddo + + do o = 1_pInt,lattice_maxNtwinFamily + index_otherFamily = sum(plastic_phenoplus_Ntwin(1:o-1_pInt,instance)) + do k = 1_pInt,plastic_phenoplus_Ntwin(o,instance) ! loop over (active) systems in other family (twin) + plastic_phenoplus_hardeningMatrix_SlipTwin(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_phenoplus_interaction_SlipTwin(lattice_interactionSlipTwin( & + sum(lattice_NslipSystem(1:f-1_pInt,phase))+j, & + sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & + phase), instance ) + enddo; enddo + + enddo; enddo + + do f = 1_pInt,lattice_maxNtwinFamily ! >>> interaction twin -- X + index_myFamily = sum(plastic_phenoplus_Ntwin(1:f-1_pInt,instance)) + do j = 1_pInt,plastic_phenoplus_Ntwin(f,instance) ! loop over (active) systems in my family (twin) + + do o = 1_pInt,lattice_maxNslipFamily + index_otherFamily = sum(plastic_phenoplus_Nslip(1:o-1_pInt,instance)) + do k = 1_pInt,plastic_phenoplus_Nslip(o,instance) ! loop over (active) systems in other family (slip) + plastic_phenoplus_hardeningMatrix_TwinSlip(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_phenoplus_interaction_TwinSlip(lattice_interactionTwinSlip( & + sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & + sum(lattice_NslipSystem(1:o-1_pInt,phase))+k, & + phase), instance ) + enddo; enddo + + do o = 1_pInt,lattice_maxNtwinFamily + index_otherFamily = sum(plastic_phenoplus_Ntwin(1:o-1_pInt,instance)) + do k = 1_pInt,plastic_phenoplus_Ntwin(o,instance) ! loop over (active) systems in other family (twin) + plastic_phenoplus_hardeningMatrix_TwinTwin(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_phenoplus_interaction_TwinTwin(lattice_interactionTwinTwin( & + sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & + sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & + phase), instance ) + enddo; enddo + + enddo; enddo + + call plastic_phenoplus_stateInit(phase,instance) + call plastic_phenoplus_aTolState(phase,instance) + endif myPhase2 + enddo initializeInstances + +end subroutine plastic_phenoplus_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the initial microstructural state for a given instance of this plasticity +!-------------------------------------------------------------------------------------------------- +subroutine plastic_phenoplus_stateInit(ph,instance) + use lattice, only: & + lattice_maxNslipFamily, & + lattice_maxNtwinFamily + use material, only: & + plasticState + + implicit none + integer(pInt), intent(in) :: & + instance, & !< number specifying the instance of the plasticity + ph + integer(pInt) :: & + i + real(pReal), dimension(plasticState(ph)%sizeState) :: & + tempState + + tempState = 0.0_pReal + do i = 1_pInt,lattice_maxNslipFamily + tempState(1+sum(plastic_phenoplus_Nslip(1:i-1,instance)) : & + sum(plastic_phenoplus_Nslip(1:i ,instance))) = & + plastic_phenoplus_tau0_slip(i,instance) + enddo + + do i = 1_pInt,lattice_maxNtwinFamily + tempState(1+sum(plastic_phenoplus_Nslip(:,instance))+& + sum(plastic_phenoplus_Ntwin(1:i-1,instance)) : & + sum(plastic_phenoplus_Nslip(:,instance))+& + sum(plastic_phenoplus_Ntwin(1:i ,instance))) = & + plastic_phenoplus_tau0_twin(i,instance) + enddo + + plasticState(ph)%state0(:,:) = spread(tempState, & ! spread single tempstate array + 2, & ! along dimension 2 + size(plasticState(ph)%state0(1,:))) ! number of copies (number of IPCs) + +end subroutine plastic_phenoplus_stateInit + + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this plasticity +!-------------------------------------------------------------------------------------------------- +subroutine plastic_phenoplus_aTolState(ph,instance) + use material, only: & + plasticState + + implicit none + integer(pInt), intent(in) :: & + instance, & !< number specifying the instance of the plasticity + ph + + plasticState(ph)%aTolState(1:plastic_phenoplus_totalNslip(instance)+ & + plastic_phenoplus_totalNtwin(instance)) = & + plastic_phenoplus_aTolResistance(instance) + plasticState(ph)%aTolState(1+plastic_phenoplus_totalNslip(instance)+ & + plastic_phenoplus_totalNtwin(instance)) = & + plastic_phenoplus_aTolShear(instance) + plasticState(ph)%aTolState(2+plastic_phenoplus_totalNslip(instance)+ & + plastic_phenoplus_totalNtwin(instance)) = & + plastic_phenoplus_aTolTwinFrac(instance) + plasticState(ph)%aTolState(3+plastic_phenoplus_totalNslip(instance)+ & + plastic_phenoplus_totalNtwin(instance): & + 2+2*(plastic_phenoplus_totalNslip(instance)+ & + plastic_phenoplus_totalNtwin(instance))) = & + plastic_phenoplus_aTolShear(instance) + +end subroutine plastic_phenoplus_aTolState + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculate push-up factors (kappa) for each voxel based on its neighbors +!-------------------------------------------------------------------------------------------------- +subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el) + use math, only: pi, & + math_mul33x33, & + math_mul3x3, & + math_transpose33, & + math_qDot, & + math_qRot, & + indeg + + use mesh, only: mesh_element, & + FE_NipNeighbors, & + FE_geomtype, & + FE_celltype, & + mesh_maxNips, & + mesh_NcpElems, & + mesh_ipNeighborhood + + use material, only: material_phase, & + material_texture, & + phase_plasticityInstance, & + mappingConstitutive, & + homogenization_maxNgrains, & + plasticState + + use lattice, only: lattice_sn, & + lattice_sd, & + lattice_qDisorientation + + !***input variables + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + orientation ! crystal orientation in quaternions + + !***local variables + integer(pInt) instance, & !my instance of this plasticity + ph, & !my phase + of, & !my spatial position in memory (offset) + textureID, & !my texture + Nneighbors, & !number of neighbors (<= 6) + n, & !neighbor index (for iterating through all neighbors) + ns, & !number of slip system + nns, & !number of slip in my neighbors + nt, & !number of twin system + me_slip, & !my slip system index + neighbor_el, & !element number of neighboring material point + neighbor_ip, & !integration point of neighboring material point + neighbor_n, & !I have no idea what is this + neighbor_of, & !spatial position in memory for this neighbor (offset) + neighbor_ph, & !neighbor's phase + ne_slip, & !slip system index for neighbor + index_kappa, & !index of pushup factors in plasticState + offset_acshear_slip, & !offset in PlasticState for the accumulative shear + j !quickly loop through slip families + + real(pReal) kappa_max, & ! + weight_sum, & !temp storage for the denominator in weighting function + mprimeavg + + + real(pReal), dimension(plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & + m_primes, & !m' between me_alpha(one) and neighbor beta(all) + me_acshear, & !temp storage for ac_shear of one particular system for me + ne_acshear !temp storage for ac_shear of one particular system for one of my neighbor + + real(pReal), dimension(3,plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & + slipNormal, & + slipDirect + + real(pReal), dimension(4) :: my_orientation, & !store my orientation + neighbor_orientation, & !store my neighbor orientation + absMisorientation + + real(pReal), dimension(FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) :: & + ne_mprimes !m' between each neighbor + + !***Get my properties + Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) + ph = mappingConstitutive(2,ipc,ip,el) !get my phase + of = mappingConstitutive(1,ipc,ip,el) !get my spatial location offset in memory + textureID = material_texture(1,ip,el) + instance = phase_plasticityInstance(ph) !get my instance based on phase ID + ns = plastic_phenoplus_totalNslip(instance) + nt = plastic_phenoplus_totalNtwin(instance) + offset_acshear_slip = ns + nt + 2_pInt + kappa_max = ns + nt + 2_pInt + ns + nt !location of kappa in plasticState + + !***gather my accumulative shear from palsticState + findMyShear: do j = 1_pInt,ns + me_acshear(j) = plasticState(ph)%state(offset_acshear_slip+j, of) + enddo findMyShear + + !***gather my orientation and slip systems + my_orientation = orientation(1:4, ipc, ip, el) + slipNormal(1:3, 1:ns) = lattice_sn(1:3, 1:ns, ph) + slipDirect(1:3, 1:ns) = lattice_sd(1:3, 1:ns, ph) + kappa_max = plastic_phenoplus_kappa_max(instance) !maximum pushups allowed (READIN) + + !***calculate kappa between me and all my neighbors + loopMySlip: do me_slip=1_pInt,ns + loopNeighbors: do n=1_pInt,Nneighbors + neighbor_el = mesh_ipNeighborhood(1,n,ip,el) + neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) + neighbor_n = 1 !It is ipc + neighbor_of = mappingConstitutive(1, neighbor_n, neighbor_ip, neighbor_el) + neighbor_ph = mappingConstitutive(2, neighbor_n, neighbor_ip, neighbor_el) + neighbor_orientation = orientation(1:4, neighbor_n, neighbor_ip, neighbor_el) !ipc is always 1. + absMisorientation = lattice_qDisorientation(my_orientation, & + neighbor_orientation, & + 0_pInt) !no need for explicit calculation of symmetry + if (ph==neighbor_ph) then !ignore voxel with different phase than phenoplus + loopNeighborSlip: do ne_slip=1_pInt,ns !assuming all have the same slip systems here + m_primes(ne_slip) = abs(math_mul3x3(slipNormal(1:3,me_slip), & + math_qRot(absMisorientation, slipNormal(1:3,ne_slip)))) & + * abs(math_mul3x3(slipDirect(1:3,me_slip), & + math_qRot(absMisorientation, slipDirect(1:3,ne_slip)))) !follow Philip's suggestions: + !only using the geometrical compatibility for now + enddo loopNeighborSlip + ne_mprimes(n) = maxval(m_primes) !pick the highest m' out of all slip combinations + else + ne_mprimes(n) = 0.0_pReal !0 for other phase + Nneighbors = Nneighbors - 1_pInt !kick neighbor with different phase out + endif + + enddo loopNeighbors + + mprimeavg = sum(ne_mprimes)/Nneighbors !average the m' from all "good" neighbors + plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal + & !calculate kappa + 2.0_pReal * & + (kappa_max - 1.0_pReal) * & + (1.0_pReal - mprimeavg) + enddo loopMySlip + +end subroutine plastic_phenoplus_microstructure + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- +subroutine plastic_phenoplus_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el) + use math, only: & + math_Plain3333to99, & + math_Mandel6to33 + use lattice, only: & + lattice_Sslip, & + lattice_Sslip_v, & + lattice_Stwin, & + lattice_Stwin_v, & + lattice_maxNslipFamily, & + lattice_maxNtwinFamily, & + lattice_NslipSystem, & + lattice_NtwinSystem, & + lattice_NnonSchmid + use material, only: & + plasticState, & + mappingConstitutive, & + phase_plasticityInstance + + implicit none + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(9,9), intent(out) :: & + dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress + + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + + integer(pInt) :: & + instance, & + nSlip, & + nTwin,index_Gamma,index_F,index_myFamily, index_kappa, & + f,i,j,k,l,m,n, & + of, & + ph + real(pReal) :: & + tau_slip_pos,tau_slip_neg, & + gdot_slip_pos,gdot_slip_neg, & + dgdot_dtauslip_pos,dgdot_dtauslip_neg, & + gdot_twin,dgdot_dtautwin,tau_twin + real(pReal), dimension(3,3,3,3) :: & + dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor + real(pReal), dimension(3,3,2) :: & + nonSchmid_tensor + + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + nSlip = plastic_phenoplus_totalNslip(instance) + nTwin = plastic_phenoplus_totalNtwin(instance) + index_Gamma = nSlip + nTwin + 1_pInt + index_F = nSlip + nTwin + 2_pInt + index_kappa = nSlip + nTwin + 2_pInt +nSlip + nTwin + + Lp = 0.0_pReal + dLp_dTstar3333 = 0.0_pReal + dLp_dTstar99 = 0.0_pReal + +!-------------------------------------------------------------------------------------------------- +! Slip part + j = 0_pInt + slipFamilies: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems: do i = 1_pInt,plastic_phenoplus_Nslip(f,instance) + j = j+1_pInt + + ! Calculation of Lp + tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_pos + nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) + nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) + do k = 1,lattice_NnonSchmid(ph) + tau_slip_pos = tau_slip_pos + plastic_phenoplus_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_neg + plastic_phenoplus_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) + nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + plastic_phenoplus_nonSchmidCoeff(k,instance)*& + lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph) + nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + plastic_phenoplus_nonSchmidCoeff(k,instance)*& + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph) + enddo + + !insert non-local effect here by modify gdot with kappa in plastic state + gdot_slip_pos = 0.5_pReal*plastic_phenoplus_gdot0_slip(instance)* & + ((abs(tau_slip_pos)/(plasticState(ph)%state(j, of)* & + plasticState(ph)%state(j+index_kappa, of))) & !in-place modification of gdot + **plastic_phenoplus_n_slip(instance))*sign(1.0_pReal,tau_slip_pos) + + gdot_slip_neg = 0.5_pReal*plastic_phenoplus_gdot0_slip(instance)* & + ((abs(tau_slip_neg)/(plasticState(ph)%state(j, of)* & + plasticState(ph)%state(j+index_kappa, of))) & !?should we make it direction aware + **plastic_phenoplus_n_slip(instance))*sign(1.0_pReal,tau_slip_neg) + + Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F + (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) + + ! Calculation of the tangent of Lp + if (abs(gdot_slip_pos) > tiny(0.0_pReal)) then + dgdot_dtauslip_pos = gdot_slip_pos*plastic_phenoplus_n_slip(instance)/tau_slip_pos + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & + dgdot_dtauslip_pos*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & + nonSchmid_tensor(m,n,1) + endif + + if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then + dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenoplus_n_slip(instance)/tau_slip_neg + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & + dgdot_dtauslip_neg*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & + nonSchmid_tensor(m,n,2) + endif + enddo slipSystems + enddo slipFamilies + +!-------------------------------------------------------------------------------------------------- +! Twinning part + j = 0_pInt + twinFamilies: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems: do i = 1_pInt,plastic_phenoplus_Ntwin(f,instance) + j = j+1_pInt + + ! Calculation of Lp + tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) + gdot_twin = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F + plastic_phenoplus_gdot0_twin(instance)*& + (abs(tau_twin)/plasticState(ph)%state(nSlip+j,of))**& + plastic_phenoplus_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin)) + Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph) + + ! Calculation of the tangent of Lp + if (abs(gdot_twin) > tiny(0.0_pReal)) then + dgdot_dtautwin = gdot_twin*plastic_phenoplus_n_twin(instance)/tau_twin + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & + dgdot_dtautwin*lattice_Stwin(k,l,index_myFamily+i,ph)* & + lattice_Stwin(m,n,index_myFamily+i,ph) + endif + enddo twinSystems + enddo twinFamilies + + dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + + +end subroutine plastic_phenoplus_LpAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +subroutine plastic_phenoplus_dotState(Tstar_v,ipc,ip,el) + use lattice, only: & + lattice_Sslip_v, & + lattice_Stwin_v, & + lattice_maxNslipFamily, & + lattice_maxNtwinFamily, & + lattice_NslipSystem, & + lattice_NtwinSystem, & + lattice_shearTwin, & + lattice_NnonSchmid + use material, only: & + material_phase, & + mappingConstitutive, & + plasticState, & + phase_plasticityInstance + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element !< microstructure state + + integer(pInt) :: & + instance,ph, & + nSlip,nTwin, & + f,i,j,k, & + index_Gamma,index_F,index_myFamily, & + offset_accshear_slip,offset_accshear_twin, & + of + real(pReal) :: & + c_SlipSlip,c_TwinSlip,c_TwinTwin, & + ssat_offset, & + tau_slip_pos,tau_slip_neg,tau_twin + + real(pReal), dimension(plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip + real(pReal), dimension(plastic_phenoplus_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + gdot_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin + + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + + nSlip = plastic_phenoplus_totalNslip(instance) + nTwin = plastic_phenoplus_totalNtwin(instance) + + index_Gamma = nSlip + nTwin + 1_pInt + index_F = nSlip + nTwin + 2_pInt + offset_accshear_slip = nSlip + nTwin + 2_pInt + offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip + plasticState(ph)%dotState(:,of) = 0.0_pReal + + +!-------------------------------------------------------------------------------------------------- +! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices + c_SlipSlip = plastic_phenoplus_h0_SlipSlip(instance)*& + (1.0_pReal + plastic_phenoplus_twinC(instance)*plasticState(ph)%state(index_F,of)**& + plastic_phenoplus_twinB(instance)) + c_TwinSlip = plastic_phenoplus_h0_TwinSlip(instance)*& + plasticState(ph)%state(index_Gamma,of)**plastic_phenoplus_twinE(instance) + c_TwinTwin = plastic_phenoplus_h0_TwinTwin(instance)*& + plasticState(ph)%state(index_F,of)**plastic_phenoplus_twinD(instance) + +!-------------------------------------------------------------------------------------------------- +! calculate left and right vectors and calculate dot gammas + ssat_offset = plastic_phenoplus_spr(instance)*sqrt(plasticState(ph)%state(index_F,of)) + j = 0_pInt + slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems1: do i = 1_pInt,plastic_phenoplus_Nslip(f,instance) + j = j+1_pInt + left_SlipSlip(j) = 1.0_pReal ! no system-dependent left part + left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part + right_SlipSlip(j) = abs(1.0_pReal-plasticState(ph)%state(j,of) / & + (plastic_phenoplus_tausat_slip(f,instance)+ssat_offset)) & + **plastic_phenoplus_a_slip(instance)& + *sign(1.0_pReal,1.0_pReal-plasticState(ph)%state(j,of) / & + (plastic_phenoplus_tausat_slip(f,instance)+ssat_offset)) + right_TwinSlip(j) = 1.0_pReal ! no system-dependent part + +!-------------------------------------------------------------------------------------------------- +! Calculation of dot gamma + tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_pos + nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph) + tau_slip_pos = tau_slip_pos + plastic_phenoplus_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k, index_myFamily+i,ph)) + tau_slip_neg = tau_slip_neg + plastic_phenoplus_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) + enddo nonSchmidSystems + gdot_slip(j) = plastic_phenoplus_gdot0_slip(instance)*0.5_pReal* & + ((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenoplus_n_slip(instance) & + +(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenoplus_n_slip(instance))& + *sign(1.0_pReal,tau_slip_pos) + enddo slipSystems1 + enddo slipFamilies1 + + + j = 0_pInt + twinFamilies1: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems1: do i = 1_pInt,plastic_phenoplus_Ntwin(f,instance) + j = j+1_pInt + left_TwinSlip(j) = 1.0_pReal ! no system-dependent left part + left_TwinTwin(j) = 1.0_pReal ! no system-dependent left part + right_SlipTwin(j) = 1.0_pReal ! no system-dependent right part + right_TwinTwin(j) = 1.0_pReal ! no system-dependent right part + +!-------------------------------------------------------------------------------------------------- +! Calculation of dot vol frac + tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) + gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F + plastic_phenoplus_gdot0_twin(instance)*& + (abs(tau_twin)/plasticState(ph)%state(nslip+j,of))**& + plastic_phenoplus_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin)) + enddo twinSystems1 + enddo twinFamilies1 + +!-------------------------------------------------------------------------------------------------- +! calculate the overall hardening based on above + j = 0_pInt + slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + slipSystems2: do i = 1_pInt,plastic_phenoplus_Nslip(f,instance) + j = j+1_pInt + plasticState(ph)%dotState(j,of) = & ! evolution of slip resistance j + c_SlipSlip * left_SlipSlip(j) * & + dot_product(plastic_phenoplus_hardeningMatrix_SlipSlip(j,1:nSlip,instance), & + right_SlipSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor + dot_product(plastic_phenoplus_hardeningMatrix_SlipTwin(j,1:nTwin,instance), & + right_SlipTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor + plasticState(ph)%dotState(index_Gamma,of) = plasticState(ph)%dotState(index_Gamma,of) + & + abs(gdot_slip(j)) + plasticState(ph)%dotState(offset_accshear_slip+j,of) = abs(gdot_slip(j)) + enddo slipSystems2 + enddo slipFamilies2 + + j = 0_pInt + twinFamilies2: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems2: do i = 1_pInt,plastic_phenoplus_Ntwin(f,instance) + j = j+1_pInt + plasticState(ph)%dotState(j+nSlip,of) = & ! evolution of twin resistance j + c_TwinSlip * left_TwinSlip(j) * & + dot_product(plastic_phenoplus_hardeningMatrix_TwinSlip(j,1:nSlip,instance), & + right_TwinSlip*abs(gdot_slip)) + & ! dot gamma_slip modulated by right-side slip factor + c_TwinTwin * left_TwinTwin(j) * & + dot_product(plastic_phenoplus_hardeningMatrix_TwinTwin(j,1:nTwin,instance), & + right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor + if (plasticState(ph)%state(index_F,of) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 + plasticState(ph)%dotState(index_F,of) = plasticState(ph)%dotState(index_F,of) + & + gdot_twin(j)/lattice_shearTwin(index_myFamily+i,ph) + plasticState(ph)%dotState(offset_accshear_twin+j,of) = abs(gdot_twin(j)) + enddo twinSystems2 + enddo twinFamilies2 + + +end subroutine plastic_phenoplus_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function plastic_phenoplus_postResults(Tstar_v,ipc,ip,el) + use material, only: & + material_phase, & + plasticState, & + mappingConstitutive, & + phase_plasticityInstance + use lattice, only: & + lattice_Sslip_v, & + lattice_Stwin_v, & + lattice_maxNslipFamily, & + lattice_maxNtwinFamily, & + lattice_NslipSystem, & + lattice_NtwinSystem, & + lattice_NnonSchmid + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element !< microstructure state + + real(pReal), dimension(plastic_phenoplus_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + plastic_phenoplus_postResults + + integer(pInt) :: & + instance,ph, of, & + nSlip,nTwin, & + o,f,i,c,j,k, & + index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily,index_kappa + real(pReal) :: & + tau_slip_pos,tau_slip_neg,tau + + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + + nSlip = plastic_phenoplus_totalNslip(instance) + nTwin = plastic_phenoplus_totalNtwin(instance) + + index_Gamma = nSlip + nTwin + 1_pInt + index_F = nSlip + nTwin + 2_pInt + index_accshear_slip = nSlip + nTwin + 2_pInt + 1_pInt + index_accshear_twin = nSlip + nTwin + 2_pInt + nSlip + 1_pInt + index_kappa = nSlip + nTwin + 2_pInt + nSlip + nTwin + 1_pInt + + plastic_phenoplus_postResults = 0.0_pReal + c = 0_pInt + + outputsLoop: do o = 1_pInt,plastic_phenoplus_Noutput(instance) + select case(plastic_phenoplus_outputID(o,instance)) + case (resistance_slip_ID) + plastic_phenoplus_postResults(c+1_pInt:c+nSlip) = plasticState(ph)%state(1:nSlip,of) + c = c + nSlip + + case (accumulatedshear_slip_ID) + plastic_phenoplus_postResults(c+1_pInt:c+nSlip) = plasticState(ph)%state(index_accshear_slip:& + index_accshear_slip+nSlip-1_pInt,of) + c = c + nSlip + + case (shearrate_slip_ID) + j = 0_pInt + slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems1: do i = 1_pInt,plastic_phenoplus_Nslip(f,instance) + j = j + 1_pInt + tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_pos + do k = 1,lattice_NnonSchmid(ph) + tau_slip_pos = tau_slip_pos + plastic_phenoplus_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_neg + plastic_phenoplus_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) + enddo + plastic_phenoplus_postResults(c+j) = plastic_phenoplus_gdot0_slip(instance)*0.5_pReal* & + ((abs(tau_slip_pos)/plasticState(ph)%state(j,of))**plastic_phenoplus_n_slip(instance) & + +(abs(tau_slip_neg)/plasticState(ph)%state(j,of))**plastic_phenoplus_n_slip(instance))& + *sign(1.0_pReal,tau_slip_pos) + + enddo slipSystems1 + enddo slipFamilies1 + c = c + nSlip + + case (resolvedstress_slip_ID) + j = 0_pInt + slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems2: do i = 1_pInt,plastic_phenoplus_Nslip(f,instance) + j = j + 1_pInt + plastic_phenoplus_postResults(c+j) = & + dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + enddo slipSystems2 + enddo slipFamilies2 + c = c + nSlip + + case (kappa_slip_ID) + plastic_phenoplus_postResults(c+1_pInt:c+nSlip) = & + plasticState(ph)%state(index_kappa:index_kappa+nSlip-1_pInt,of) + c = c + nSlip + + case (totalshear_ID) + plastic_phenoplus_postResults(c+1_pInt) = & + plasticState(ph)%state(index_Gamma,of) + c = c + 1_pInt + + case (resistance_twin_ID) + plastic_phenoplus_postResults(c+1_pInt:c+nTwin) = & + plasticState(ph)%state(1_pInt+nSlip:nTwin+nSlip-1_pInt,of) + c = c + nTwin + + case (accumulatedshear_twin_ID) + plastic_phenoplus_postResults(c+1_pInt:c+nTwin) = & + plasticState(ph)%state(index_accshear_twin:index_accshear_twin+nTwin-1_pInt,of) + c = c + nTwin + + case (shearrate_twin_ID) + j = 0_pInt + twinFamilies1: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems1: do i = 1_pInt,plastic_phenoplus_Ntwin(f,instance) + j = j + 1_pInt + tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) + plastic_phenoplus_postResults(c+j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F + plastic_phenoplus_gdot0_twin(instance)*& + (abs(tau)/plasticState(ph)%state(j+nSlip,of))**& + plastic_phenoplus_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau)) + enddo twinSystems1 + enddo twinFamilies1 + c = c + nTwin + + case (resolvedstress_twin_ID) + j = 0_pInt + twinFamilies2: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems2: do i = 1_pInt,plastic_phenoplus_Ntwin(f,instance) + j = j + 1_pInt + plastic_phenoplus_postResults(c+j) = & + dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) + enddo twinSystems2 + enddo twinFamilies2 + c = c + nTwin + + case (totalvolfrac_twin_ID) + plastic_phenoplus_postResults(c+1_pInt) = plasticState(ph)%state(index_F,of) + c = c + 1_pInt + + end select + enddo outputsLoop + +end function plastic_phenoplus_postResults + +end module plastic_phenoplus diff --git a/code/plastic_phenopowerlaw.f90 b/code/plastic_phenopowerlaw.f90 index 8c585630e..220e11907 100644 --- a/code/plastic_phenopowerlaw.f90 +++ b/code/plastic_phenopowerlaw.f90 @@ -3,7 +3,7 @@ !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw +!> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw !! fitting !-------------------------------------------------------------------------------------------------- module plastic_phenopowerlaw @@ -19,11 +19,11 @@ module plastic_phenopowerlaw integer(pInt), dimension(:,:), allocatable, target, public :: & plastic_phenopowerlaw_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & + character(len=64), dimension(:,:), allocatable, target, public :: & plastic_phenopowerlaw_output !< name of each post result output integer(pInt), dimension(:), allocatable, target, public :: & - plastic_phenopowerlaw_Noutput !< number of outputs per instance of this constitution + plastic_phenopowerlaw_Noutput !< number of outputs per instance of this constitution integer(pInt), dimension(:), allocatable, public, protected :: & plastic_phenopowerlaw_totalNslip, & !< no. of slip system used in simulation @@ -53,9 +53,9 @@ module plastic_phenopowerlaw plastic_phenopowerlaw_aTolShear, & plastic_phenopowerlaw_aTolTwinfrac, & plastic_phenopowerlaw_aTolTransfrac, & - plastic_phenopowerlaw_Cnuc, & !< coefficient for strain-induced martensite nucleation - plastic_phenopowerlaw_Cdwp, & !< coefficient for double well potential - plastic_phenopowerlaw_Cgro, & !< coefficient for stress-assisted martensite growth + plastic_phenopowerlaw_Cnuc, & !< coefficient for strain-induced martensite nucleation + plastic_phenopowerlaw_Cdwp, & !< coefficient for double well potential + plastic_phenopowerlaw_Cgro, & !< coefficient for stress-assisted martensite growth plastic_phenopowerlaw_deltaG !< free energy difference between austensite and martensite [MPa] real(pReal), dimension(:,:), allocatable, private :: & @@ -75,7 +75,7 @@ module plastic_phenopowerlaw plastic_phenopowerlaw_hardeningMatrix_TwinSlip, & plastic_phenopowerlaw_hardeningMatrix_TwinTwin - enum, bind(c) + enum, bind(c) enumerator :: undefined_ID, & resistance_slip_ID, & accumulatedshear_slip_ID, & @@ -88,9 +88,9 @@ module plastic_phenopowerlaw resolvedstress_twin_ID, & totalvolfrac_twin_ID end enum - integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & plastic_phenopowerlaw_outputID !< ID of each post result output - + public :: & plastic_phenopowerlaw_init, & plastic_phenopowerlaw_LpAndItsTangent, & @@ -163,14 +163,14 @@ subroutine plastic_phenopowerlaw_init(fileUnit) tag = '', & line = '' real(pReal), dimension(:), allocatable :: tempPerSlip - - mainProcess: if (worldrank == 0) then + + mainProcess: if (worldrank == 0) then write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" endif mainProcess - + maxNinstance = int(count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID),pInt) if (maxNinstance == 0_pInt) return @@ -239,7 +239,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) if (IO_isBlank(line)) cycle ! skip empty lines if (IO_getTag(line,'<','>') /= '') then ! stop at next part line = IO_read(fileUnit, .true.) ! reset IO_read - exit + exit endif if (IO_getTag(line,'[',']') /= '') then ! next phase phase = phase + 1_pInt ! advance phase section counter @@ -451,7 +451,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) sanityChecks: do phase = 1_pInt, size(phase_plasticity) myPhase: if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then - instance = phase_plasticityInstance(phase) + instance = phase_plasticityInstance(phase) plastic_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance) = & min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested plastic_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance)) @@ -537,12 +537,12 @@ subroutine plastic_phenopowerlaw_init(fileUnit) mySize = 1_pInt case default end select - + outputFound: if (mySize > 0_pInt) then plastic_phenopowerlaw_sizePostResult(o,instance) = mySize plastic_phenopowerlaw_sizePostResults(instance) = plastic_phenopowerlaw_sizePostResults(instance) + mySize endif outputFound - enddo outputsLoop + enddo outputsLoop !-------------------------------------------------------------------------------------------------- ! allocate state arrays sizeState = plastic_phenopowerlaw_totalNslip(instance) & ! s_slip @@ -550,7 +550,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) + 2_pInt & ! sum(gamma) + sum(f) + plastic_phenopowerlaw_totalNslip(instance) & ! accshear_slip + plastic_phenopowerlaw_totalNtwin(instance) ! accshear_twin - + sizeDotState = sizeState sizeDeltaState = 0_pInt plasticState(phase)%sizeState = sizeState @@ -579,13 +579,13 @@ subroutine plastic_phenopowerlaw_init(fileUnit) allocate(plasticState(phase)%RK4dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) if (any(numerics_integrator == 5_pInt)) & allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase), source=0.0_pReal) - + offset_slip = plasticState(phase)%nSlip+plasticState(phase)%nTwin+2_pInt plasticState(phase)%slipRate => & plasticState(phase)%dotState(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase) plasticState(phase)%accumulatedSlip => & plasticState(phase)%state(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NipcMyPhase) - + do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X index_myFamily = sum(plastic_phenopowerlaw_Nslip(1:f-1_pInt,instance)) do j = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) ! loop over (active) systems in my family (slip) @@ -598,7 +598,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) sum(lattice_NslipSystem(1:o-1,phase))+k, & phase), instance ) enddo; enddo - + do o = 1_pInt,lattice_maxNtwinFamily index_otherFamily = sum(plastic_phenopowerlaw_Ntwin(1:o-1_pInt,instance)) do k = 1_pInt,plastic_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin) @@ -608,13 +608,13 @@ subroutine plastic_phenopowerlaw_init(fileUnit) sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & phase), instance ) enddo; enddo - + enddo; enddo - + do f = 1_pInt,lattice_maxNtwinFamily ! >>> interaction twin -- X index_myFamily = sum(plastic_phenopowerlaw_Ntwin(1:f-1_pInt,instance)) do j = 1_pInt,plastic_phenopowerlaw_Ntwin(f,instance) ! loop over (active) systems in my family (twin) - + do o = 1_pInt,lattice_maxNslipFamily index_otherFamily = sum(plastic_phenopowerlaw_Nslip(1:o-1_pInt,instance)) do k = 1_pInt,plastic_phenopowerlaw_Nslip(o,instance) ! loop over (active) systems in other family (slip) @@ -624,7 +624,7 @@ subroutine plastic_phenopowerlaw_init(fileUnit) sum(lattice_NslipSystem(1:o-1_pInt,phase))+k, & phase), instance ) enddo; enddo - + do o = 1_pInt,lattice_maxNtwinFamily index_otherFamily = sum(plastic_phenopowerlaw_Ntwin(1:o-1_pInt,instance)) do k = 1_pInt,plastic_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin) @@ -634,12 +634,12 @@ subroutine plastic_phenopowerlaw_init(fileUnit) sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & phase), instance ) enddo; enddo - + enddo; enddo call plastic_phenopowerlaw_stateInit(phase,instance) call plastic_phenopowerlaw_aTolState(phase,instance) - endif myPhase2 + endif myPhase2 enddo initializeInstances end subroutine plastic_phenopowerlaw_init @@ -654,11 +654,11 @@ subroutine plastic_phenopowerlaw_stateInit(ph,instance) lattice_maxNtwinFamily use material, only: & plasticState - + implicit none integer(pInt), intent(in) :: & instance, & !< number specifying the instance of the plasticity - ph + ph integer(pInt) :: & i real(pReal), dimension(plasticState(ph)%sizeState) :: & @@ -694,7 +694,7 @@ subroutine plastic_phenopowerlaw_aTolState(ph,instance) plasticState implicit none - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & instance, & !< number specifying the instance of the plasticity ph @@ -752,7 +752,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation integer(pInt) :: & - instance, & + instance, & nSlip, & nTwin,index_Gamma,index_F,index_myFamily, & f,i,j,k,l,m,n, & @@ -767,7 +767,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor real(pReal), dimension(3,3,2) :: & nonSchmid_tensor - + of = mappingConstitutive(1,ipc,ip,el) ph = mappingConstitutive(2,ipc,ip,el) instance = phase_plasticityInstance(ph) @@ -787,13 +787,13 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family slipSystems: do i = 1_pInt,plastic_phenopowerlaw_Nslip(f,instance) j = j+1_pInt - + ! Calculation of Lp tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) tau_slip_neg = tau_slip_pos nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) - do k = 1,lattice_NnonSchmid(ph) + do k = 1,lattice_NnonSchmid(ph) tau_slip_pos = tau_slip_pos + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) tau_slip_neg = tau_slip_neg + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & @@ -810,7 +810,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, gdot_slip_neg = 0.5_pReal*plastic_phenopowerlaw_gdot0_slip(instance)* & ((abs(tau_slip_neg)/(plasticState(ph)%state(j,of))) & **plastic_phenopowerlaw_n_slip(instance))*sign(1.0_pReal,tau_slip_neg) - + Lp = Lp + (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F (gdot_slip_pos+gdot_slip_neg)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) @@ -822,7 +822,7 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, dgdot_dtauslip_pos*lattice_Sslip(k,l,1,index_myFamily+i,ph)* & nonSchmid_tensor(m,n,1) endif - + if (abs(gdot_slip_neg) > tiny(0.0_pReal)) then dgdot_dtauslip_neg = gdot_slip_neg*plastic_phenopowerlaw_n_slip(instance)/tau_slip_neg forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & @@ -842,11 +842,11 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip, j = j+1_pInt ! Calculation of Lp - tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) + tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) gdot_twin = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F plastic_phenopowerlaw_gdot0_twin(instance)*& (abs(tau_twin)/plasticState(ph)%state(nSlip+j,of))**& - plastic_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin)) + plastic_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin)) Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph) ! Calculation of the tangent of Lp @@ -877,13 +877,13 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) lattice_NslipSystem, & lattice_NtwinSystem, & lattice_shearTwin, & - lattice_NnonSchmid + lattice_NnonSchmid use material, only: & material_phase, & mappingConstitutive, & plasticState, & phase_plasticityInstance - + implicit none real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation @@ -908,11 +908,11 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) gdot_slip,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip real(pReal), dimension(plastic_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & gdot_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin - + of = mappingConstitutive(1,ipc,ip,el) ph = mappingConstitutive(2,ipc,ip,el) instance = phase_plasticityInstance(ph) - + nSlip = plastic_phenopowerlaw_totalNslip(instance) nTwin = plastic_phenopowerlaw_totalNtwin(instance) @@ -921,7 +921,7 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) offset_accshear_slip = nSlip + nTwin + 2_pInt offset_accshear_twin = nSlip + nTwin + 2_pInt + nSlip plasticState(ph)%dotState(:,of) = 0.0_pReal - + !-------------------------------------------------------------------------------------------------- ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices @@ -949,12 +949,12 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) *sign(1.0_pReal,1.0_pReal-plasticState(ph)%state(j,of) / & (plastic_phenopowerlaw_tausat_slip(f,instance)+ssat_offset)) right_TwinSlip(j) = 1.0_pReal ! no system-dependent part - + !-------------------------------------------------------------------------------------------------- -! Calculation of dot gamma +! Calculation of dot gamma tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) tau_slip_neg = tau_slip_pos - nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph) + nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph) tau_slip_pos = tau_slip_pos + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k, index_myFamily+i,ph)) tau_slip_neg = tau_slip_neg + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & @@ -963,7 +963,7 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) gdot_slip(j) = plastic_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* & ((abs(tau_slip_pos)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance) & +(abs(tau_slip_neg)/(plasticState(ph)%state(j,of)))**plastic_phenopowerlaw_n_slip(instance))& - *sign(1.0_pReal,tau_slip_pos) + *sign(1.0_pReal,tau_slip_pos) enddo slipSystems1 enddo slipFamilies1 @@ -980,7 +980,7 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) !-------------------------------------------------------------------------------------------------- ! Calculation of dot vol frac - tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) + tau_twin = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,ph)) gdot_twin(j) = (1.0_pReal-plasticState(ph)%state(index_F,of))*& ! 1-F plastic_phenopowerlaw_gdot0_twin(instance)*& (abs(tau_twin)/plasticState(ph)%state(nslip+j,of))**& @@ -1025,7 +1025,7 @@ subroutine plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) enddo twinSystems2 enddo twinFamilies2 - + end subroutine plastic_phenopowerlaw_dotState !-------------------------------------------------------------------------------------------------- @@ -1044,7 +1044,7 @@ function plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) lattice_maxNtwinFamily, & lattice_NslipSystem, & lattice_NtwinSystem, & - lattice_NnonSchmid + lattice_NnonSchmid implicit none real(pReal), dimension(6), intent(in) :: & @@ -1061,7 +1061,7 @@ function plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) instance,ph, of, & nSlip,nTwin, & o,f,i,c,j,k, & - index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily + index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily real(pReal) :: & tau_slip_pos,tau_slip_neg,tau @@ -1099,7 +1099,7 @@ function plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) j = j + 1_pInt tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) tau_slip_neg = tau_slip_pos - do k = 1,lattice_NnonSchmid(ph) + do k = 1,lattice_NnonSchmid(ph) tau_slip_pos = tau_slip_pos + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) tau_slip_neg = tau_slip_neg + plastic_phenopowerlaw_nonSchmidCoeff(k,instance)* &