From dd23bec9aa9c1a3c32f2907c728a645bf1038858 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Dec 2020 11:03:13 +0100 Subject: [PATCH] avoid global variables --- src/constitutive.f90 | 9 ++++-- src/constitutive_mech.f90 | 62 +++++++++++++++---------------------- src/homogenization_mech.f90 | 24 ++++++-------- 3 files changed, 41 insertions(+), 54 deletions(-) diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 29fca2b33..667a23127 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -44,8 +44,6 @@ module constitutive type(rotation), dimension(:,:,:), allocatable :: & crystallite_orientation !< current orientation - real(pReal), dimension(:,:,:,:,:), allocatable, public :: & - crystallite_P !< 1st Piola-Kirchhoff stress per grain type :: tTensorContainer real(pReal), dimension(:,:,:), allocatable :: data @@ -194,6 +192,11 @@ module constitutive real(pReal), dimension(3,3) :: F_e end function constitutive_mech_getF_e + module function constitutive_mech_getP(co,ip,el) result(P) + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: P + end function constitutive_mech_getP + module function constitutive_thermal_T(co,ip,el) result(T) integer, intent(in) :: co, ip, el real(pReal) :: T @@ -411,6 +414,7 @@ module constitutive constitutive_restartRead, & integrateSourceState, & constitutive_mech_setF, & + constitutive_mech_getP, & constitutive_mech_getLp, & constitutive_mech_getF, & constitutive_mech_getS, & @@ -877,7 +881,6 @@ subroutine crystallite_init iMax = discretization_nIPs eMax = discretization_Nelems - allocate(crystallite_P(3,3,cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_orientation(cMax,iMax,eMax)) num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index a1256a15d..f5a5cd0a2 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -24,6 +24,7 @@ submodule(constitutive) constitutive_mech constitutive_mech_Li, & constitutive_mech_Lp, & constitutive_mech_S, & + constitutive_mech_P, & ! converged value at end of last solver increment constitutive_mech_Fi0, & constitutive_mech_Fp0, & @@ -363,6 +364,7 @@ module subroutine mech_init allocate(constitutive_mech_Lp0(phases%length)) allocate(constitutive_mech_Lp(phases%length)) allocate(constitutive_mech_S(phases%length)) + allocate(constitutive_mech_P(phases%length)) allocate(constitutive_mech_S0(phases%length)) allocate(constitutive_mech_partitionedS0(phases%length)) @@ -383,6 +385,7 @@ module subroutine mech_init allocate(constitutive_mech_Lp0(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_Lp(ph)%data(3,3,Nconstituents)) allocate(constitutive_mech_S(ph)%data(3,3,Nconstituents),source=0.0_pReal) + allocate(constitutive_mech_P(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_S0(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_partitionedS0(ph)%data(3,3,Nconstituents),source=0.0_pReal) allocate(constitutive_mech_F(ph)%data(3,3,Nconstituents)) @@ -1027,7 +1030,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken) call math_invert33(Fp_new,devNull,error,invFp_new) if (error) return ! error - crystallite_P (1:3,1:3,co,ip,el) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) + constitutive_mech_P(ph)%data(1:3,1:3,me) = matmul(matmul(F,invFp_new),matmul(S,transpose(invFp_new))) constitutive_mech_S(ph)%data(1:3,1:3,me) = S constitutive_mech_Lp(ph)%data(1:3,1:3,me) = Lpguess constitutive_mech_Li(ph)%data(1:3,1:3,me) = Liguess @@ -1381,29 +1384,28 @@ subroutine crystallite_results(group,ph) select case (output_constituent(ph)%label(ou)) case('F') - call results_writeDataset(group//'/mechanics/',constitutive_mech_F(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_F(ph)%data,'F',& 'deformation gradient','1') case('F_e') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Fe(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fe(ph)%data,'F_e',& 'elastic deformation gradient','1') case('F_p') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Fp(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fp(ph)%data,'F_p', & 'plastic deformation gradient','1') case('F_i') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Fi(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Fi(ph)%data,'F_i', & 'inelastic deformation gradient','1') case('L_p') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Lp(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Lp(ph)%data,'L_p', & 'plastic velocity gradient','1/s') case('L_i') - call results_writeDataset(group//'/mechanics/',constitutive_mech_Li(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_Li(ph)%data,'L_i', & 'inelastic velocity gradient','1/s') case('P') - selected_tensors = select_tensors(crystallite_P,ph) - call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_P(ph)%data,'P', & 'First Piola-Kirchhoff stress','Pa') case('S') - call results_writeDataset(group//'/mechanics/',constitutive_mech_S(ph)%data,output_constituent(ph)%label(ou),& + call results_writeDataset(group//'/mechanics/',constitutive_mech_S(ph)%data,'S', & 'Second Piola-Kirchhoff stress','Pa') case('O') select case(lattice_structure(ph)) @@ -1430,33 +1432,6 @@ subroutine crystallite_results(group,ph) contains - !------------------------------------------------------------------------------------------------ - !> @brief select tensors for output - !------------------------------------------------------------------------------------------------ - function select_tensors(dataset,ph) - - integer, intent(in) :: ph - real(pReal), dimension(:,:,:,:,:), intent(in) :: dataset - real(pReal), allocatable, dimension(:,:,:) :: select_tensors - integer :: el,ip,co,j - - allocate(select_tensors(3,3,count(material_phaseAt==ph)*discretization_nIPs)) - - j=0 - do el = 1, size(material_phaseAt,2) - do ip = 1, discretization_nIPs - do co = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains - if (material_phaseAt(co,el) == ph) then - j = j + 1 - select_tensors(1:3,1:3,j) = dataset(1:3,1:3,co,ip,el) - endif - enddo - enddo - enddo - - end function select_tensors - - !-------------------------------------------------------------------------------------------------- !> @brief select rotations for output !-------------------------------------------------------------------------------------------------- @@ -1918,6 +1893,19 @@ module function constitutive_mech_getF_e(co,ip,el) result(F_e) end function constitutive_mech_getF_e + +! getter for non-mech (e.g. thermal) +module function constitutive_mech_getP(co,ip,el) result(P) + + integer, intent(in) :: co, ip, el + real(pReal), dimension(3,3) :: P + + + P = constitutive_mech_P(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + +end function constitutive_mech_getP + + ! setter for homogenization module subroutine constitutive_mech_setF(F,co,ip,el) diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 1d0942f3e..783d08dd1 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -111,7 +111,7 @@ module subroutine mech_partition(subF,ip,el) integer, intent(in) :: & ip, & !< integration point el !< element number - + integer :: co real(pReal) :: F(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) @@ -149,35 +149,36 @@ module subroutine mech_homogenize(dt,ip,el) integer :: co,ce real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) ce = (el-1)* discretization_nIPs + ip chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el) + homogenization_P(1:3,1:3,ce) = constitutive_mech_getP(1,ip,el) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,1,ip,el) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) + Ps(:,:,co) = constitutive_mech_getP(co,ip,el) enddo call mech_isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& - crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - dPdFs, & + Ps,dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) case (HOMOGENIZATION_RGC_ID) chosenHomogenization do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el) + Ps(:,:,co) = constitutive_mech_getP(co,ip,el) enddo call mech_RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& - crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - dPdFs, & + Ps,dPdFs, & homogenization_typeInstance(material_homogenizationAt(el))) end select chosenHomogenization @@ -203,21 +204,16 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy) integer :: co real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el))) real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el) Fs(:,:,co) = constitutive_mech_getF(co,ip,el) + Ps(:,:,co) = constitutive_mech_getP(co,ip,el) enddo - doneAndHappy = & - mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), & - Fs, & - subF,& - subdt, & - dPdFs, & - ip, & - el) + doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el) else doneAndHappy = .true. endif