From a00c6743c3c65bb350996d7622ec198330c01fec Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 26 Oct 2021 11:48:54 +0200 Subject: [PATCH 01/28] symbolic notation in numerics.yaml - p_i: integration order - p_s: shape function order ensure working combination (p_s = p_i: full integration, p_s = p_i+1: reduced integration) --- PRIVATE | 2 +- src/IO.f90 | 7 ++++++- src/mesh/FEM_utilities.f90 | 20 +++++++++++++++----- src/mesh/discretization_mesh.f90 | 8 ++++---- src/mesh/mesh_mech_FEM.f90 | 12 ++++++------ 5 files changed, 32 insertions(+), 17 deletions(-) diff --git a/PRIVATE b/PRIVATE index fabe69749..0ff1d7540 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit fabe69749425e8a7aceb3b7c2758b40d97d8b809 +Subproject commit 0ff1d7540fede9611d46d2885bebbbe60dcbbfb0 diff --git a/src/IO.f90 b/src/IO.f90 index cd7c09c75..717493006 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -432,7 +432,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'Nconstituents mismatch between homogenization and material' !-------------------------------------------------------------------------------------------------- -! material error messages and related messages in mesh +! material error messages and related messages in geometry case (150) msg = 'index out of bounds' case (153) @@ -499,6 +499,11 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (710) msg = 'Closing quotation mark missing in string' +!------------------------------------------------------------------------------------------------- +! errors related to the mesh solver + case (821) + msg = 'order not supported' + !------------------------------------------------------------------------------------------------- ! errors related to the grid solver case (831) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 64bcc3896..a3856ccaa 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -19,6 +19,7 @@ module FEM_utilities use IO use discretization_mesh use homogenization + use FEM_quadrature implicit none private @@ -29,8 +30,8 @@ module FEM_utilities !-------------------------------------------------------------------------------------------------- ! field labels information - character(len=*), parameter, public :: & - FIELD_MECH_label = 'mechanical' + character(len=*), parameter, public :: & + FIELD_MECH_label = 'mechanical' enum, bind(c); enumerator :: & FIELD_UNDEFINED_ID, & @@ -86,7 +87,9 @@ subroutine FEM_utilities_init class(tNode), pointer :: & num_mesh, & debug_mesh ! pointer to mesh debug options - integer :: structOrder !< order of displacement shape functions + integer :: & + p_s, & !< order of shape functions + p_i !< integration order (quadrature rule) character(len=*), parameter :: & PETSCDEBUG = ' -snes_view -snes_monitor ' PetscErrorCode :: ierr @@ -96,7 +99,14 @@ subroutine FEM_utilities_init print'(/,a)', ' <<<+- FEM_utilities init -+>>>' num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - structOrder = num_mesh%get_asInt('structOrder', defaultVal = 2) + + p_s = num_mesh%get_asInt('p_s',defaultVal = 2) + p_i = num_mesh%get_asInt('p_i',defaultVal = p_s) + + if (p_s < 1_pInt .or. p_s > size(FEM_nQuadrature,2)) & + call IO_error(821,ext_msg='shape function order (p_s) out of bounds') + if (p_i < max(1_pInt,p_s-1_pInt) .or. p_i > p_s) & + call IO_error(821,ext_msg='integration order (p_i) out of bounds') debug_mesh => config_debug%get('mesh',defaultVal=emptyList) debugPETSc = debug_mesh%contains('PETSc') @@ -119,7 +129,7 @@ subroutine FEM_utilities_init CHKERRQ(ierr) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asString('PETSc_options',defaultVal=''),ierr) CHKERRQ(ierr) - write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder + write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr) CHKERRQ(ierr) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index e5f989484..f1d38760d 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -85,7 +85,7 @@ subroutine discretization_mesh_init(restart) materialAt class(tNode), pointer :: & num_mesh - integer :: integrationOrder !< order of quadrature rule required + integer :: p_i !< integration order (quadrature rule) type(tvec) :: coords_node0 print'(/,a)', ' <<<+- discretization_mesh init -+>>>' @@ -93,7 +93,7 @@ subroutine discretization_mesh_init(restart) !-------------------------------------------------------------------------------- ! read numerics parameter num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + p_i = num_mesh%get_asInt('p_i',defaultVal = 2) !--------------------------------------------------------------------------------- ! read debug parameters @@ -150,9 +150,9 @@ subroutine discretization_mesh_init(restart) call VecGetArrayF90(coords_node0, mesh_node0_temp,ierr) CHKERRQ(ierr) - mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder) + mesh_maxNips = FEM_nQuadrature(dimPlex,p_i) - call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) + call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,p_i)%p) call mesh_FEM_build_ipVolumes(dimPlex) allocate(materialAt(mesh_NcpElems)) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index d6d314a42..58429992c 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -41,7 +41,7 @@ module mesh_mechanical_FEM type, private :: tNumerics integer :: & - integrationOrder, & !< order of quadrature rule required + p_i, & !< integration order itmax logical :: & BBarStabilisation @@ -118,7 +118,7 @@ subroutine FEM_mechanical_init(fieldBC) !----------------------------------------------------------------------------- ! read numerical parametes and do sanity checks num_mesh => config_numerics%get('mesh',defaultVal=emptyDict) - num%integrationOrder = num_mesh%get_asInt('integrationorder',defaultVal = 2) + num%p_i = num_mesh%get_asInt('p_i',defaultVal = 2) num%itmax = num_mesh%get_asInt('itmax',defaultVal=250) num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.) num%eps_struct_atol = num_mesh%get_asFloat('eps_struct_atol', defaultVal = 1.0e-10_pReal) @@ -135,9 +135,9 @@ subroutine FEM_mechanical_init(fieldBC) !-------------------------------------------------------------------------------------------------- ! Setup FEM mech discretization - qPoints = FEM_quadrature_points( dimPlex,num%integrationOrder)%p - qWeights = FEM_quadrature_weights(dimPlex,num%integrationOrder)%p - nQuadrature = FEM_nQuadrature( dimPlex,num%integrationOrder) + qPoints = FEM_quadrature_points( dimPlex,num%p_i)%p + qWeights = FEM_quadrature_weights(dimPlex,num%p_i)%p + nQuadrature = FEM_nQuadrature( dimPlex,num%p_i) qPointsP => qPoints qWeightsP => qWeights call PetscQuadratureCreate(PETSC_COMM_SELF,mechQuad,ierr); CHKERRQ(ierr) @@ -146,7 +146,7 @@ subroutine FEM_mechanical_init(fieldBC) call PetscQuadratureSetData(mechQuad,dimPlex,nc,nQuadrature,qPointsP,qWeightsP,ierr) CHKERRQ(ierr) call PetscFECreateDefault(PETSC_COMM_SELF,dimPlex,nc,PETSC_TRUE,prefix, & - num%integrationOrder,mechFE,ierr); CHKERRQ(ierr) + num%p_i,mechFE,ierr); CHKERRQ(ierr) call PetscFESetQuadrature(mechFE,mechQuad,ierr); CHKERRQ(ierr) call PetscFEGetDimension(mechFE,nBasis,ierr); CHKERRQ(ierr) nBasis = nBasis/nc From b4ff89f3040e3a96d9616e0a0e20cfef1662bb8b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 26 Oct 2021 13:06:55 +0200 Subject: [PATCH 02/28] improved test --- PRIVATE | 2 +- src/mesh/mesh_mech_FEM.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/PRIVATE b/PRIVATE index 0ff1d7540..bb73a35f3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 0ff1d7540fede9611d46d2885bebbbe60dcbbfb0 +Subproject commit bb73a35f334831eed10d0e293aafea6c0900b7bd diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 58429992c..496a82fc5 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -41,7 +41,7 @@ module mesh_mechanical_FEM type, private :: tNumerics integer :: & - p_i, & !< integration order + p_i, & !< integration order (quadrature rule) itmax logical :: & BBarStabilisation From f0a33b452fa50b281a0c98fe1300fe0f3303770f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 28 Oct 2021 07:05:14 +0200 Subject: [PATCH 03/28] unifying style --- src/prec.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/prec.f90 b/src/prec.f90 index 2e67ae76a..d6d161a94 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -101,7 +101,7 @@ logical elemental pure function dEq(a,b,tol) dEq = abs(a-b) <= tol else dEq = abs(a-b) <= PREAL_EPSILON * maxval(abs([a,b])) - endif + end if end function dEq @@ -139,7 +139,7 @@ logical elemental pure function dEq0(a,tol) dEq0 = abs(a) <= tol else dEq0 = abs(a) <= PREAL_MIN * 10.0_pReal - endif + end if end function dEq0 @@ -178,7 +178,7 @@ logical elemental pure function cEq(a,b,tol) cEq = abs(a-b) <= tol else cEq = abs(a-b) <= PREAL_EPSILON * maxval(abs([a,b])) - endif + end if end function cEq From bb757cf82e27f1bd72c1d5a243e71d58a55fad98 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 30 Oct 2021 10:36:47 +0200 Subject: [PATCH 04/28] style adjustments --- ...phase_mechanical_plastic_dislotungsten.f90 | 191 +++++++++--------- 1 file changed, 97 insertions(+), 94 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 1e5e120e5..fce013bb2 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -148,7 +148,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) else prm%P_nS_pos = prm%P_sl prm%P_nS_neg = prm%P_sl - endif + end if prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), & phase_lattice(ph)) @@ -209,13 +209,13 @@ module function plastic_dislotungsten_init() result(myPlasticity) if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' else slipActive - rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray + rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, & prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & source = emptyRealArray) allocate(prm%forestProjection(0,0)) allocate(prm%h_sl_sl (0,0)) - endif slipActive + end if slipActive !-------------------------------------------------------------------------------------------------- ! allocate state arrays @@ -258,7 +258,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) ! exit if any parameter is out of range if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(dislotungsten)') - enddo + end do end function plastic_dislotungsten_init @@ -267,7 +267,7 @@ end function plastic_dislotungsten_init !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,ph,en) + Mp,T,ph,en) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -287,19 +287,20 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, & dot_gamma_pos,dot_gamma_neg, & ddot_gamma_dtau_pos,ddot_gamma_dtau_neg + Lp = 0.0_pReal dLp_dMp = 0.0_pReal associate(prm => param(ph)) - call kinetics(Mp,T,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) - do i = 1, prm%sum_N_sl - Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) - forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & - + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & - + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) - enddo + call kinetics(Mp,T,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg) + do i = 1, prm%sum_N_sl + Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i) + forall (k=1:3,l=1:3,m=1:3,n=1:3) & + dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) & + + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i) + end do end associate @@ -328,35 +329,36 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot_rho_dip_climb, & d_hat + associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph)) - call kinetics(Mp,T,ph,en,& - dot_gamma_pos,dot_gamma_neg, & - tau_pos_out = tau_pos,tau_neg_out = tau_neg) + call kinetics(Mp,T,ph,en,& + dot_gamma_pos,dot_gamma_neg, & + tau_pos_out = tau_pos,tau_neg_out = tau_neg) - dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) + dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) - where(dEq0(tau_pos)) ! ToDo: use avg of +/- - dot_rho_dip_formation = 0.0_pReal - dot_rho_dip_climb = 0.0_pReal - else where - d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & ! ToDo: use avg of +/- - prm%d_caron, & ! lower limit - dst%Lambda_sl(:,en)) ! upper limit - dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & - 0.0_pReal, & - prm%dipoleformation) - v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & - * (1.0_pReal/(d_hat+prm%d_caron)) - dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? - end where + where(dEq0(tau_pos)) ! ToDo: use avg of +/- + dot_rho_dip_formation = 0.0_pReal + dot_rho_dip_climb = 0.0_pReal + else where + d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & ! ToDo: use avg of +/- + prm%d_caron, & ! lower limit + dst%Lambda_sl(:,en)) ! upper limit + dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & + 0.0_pReal, & + prm%dipoleformation) + v_cl = (3.0_pReal*prm%mu*prm%D_0*exp(-prm%Q_cl/(kB*T))*prm%f_at/(2.0_pReal*PI*kB*T)) & + * (1.0_pReal/(d_hat+prm%d_caron)) + dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(d_hat-prm%d_caron) ! ToDo: Discuss with Franz: Stress dependency? + end where - dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication - - dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges - dot%rho_dip(:,en) = dot_rho_dip_formation & - - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole - - dot_rho_dip_climb + dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication + - dot_rho_dip_formation & + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_mob(:,en)*dot%gamma_sl(:,en) ! Spontaneous annihilation of 2 edges + dot%rho_dip(:,en) = dot_rho_dip_formation & + - (2.0_pReal*prm%d_caron)/prm%b_sl*stt%rho_dip(:,en)*dot%gamma_sl(:,en) & ! Spontaneous annihilation of an edge with a dipole + - dot_rho_dip_climb end associate @@ -368,7 +370,7 @@ end subroutine dislotungsten_dotState !-------------------------------------------------------------------------------------------------- module subroutine dislotungsten_dependentState(ph,en) - integer, intent(in) :: & + integer, intent(in) :: & ph, & en @@ -423,7 +425,7 @@ module subroutine plastic_dislotungsten_results(ph,group) 'threshold stress for slip','Pa',prm%systems_sl) end select - enddo + end do end associate @@ -456,6 +458,7 @@ pure subroutine kinetics(Mp,T,ph,en, & ddot_gamma_dtau_neg, & tau_pos_out, & tau_neg_out + real(pReal), dimension(param(ph)%sum_N_sl) :: & StressRatio, & StressRatio_p,StressRatio_pminus1, & @@ -464,80 +467,80 @@ pure subroutine kinetics(Mp,T,ph,en, & t_n, t_k, dtk,dtn integer :: j + associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - do j = 1, prm%sum_N_sl - tau_pos(j) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,j)) - tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j)) - enddo + do j = 1, prm%sum_N_sl + tau_pos(j) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,j)) + tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j)) + end do + if (present(tau_pos_out)) tau_pos_out = tau_pos + if (present(tau_neg_out)) tau_neg_out = tau_neg - if (present(tau_pos_out)) tau_pos_out = tau_pos - if (present(tau_neg_out)) tau_neg_out = tau_neg + associate(BoltzmannRatio => prm%Q_s/(kB*T), & + dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & + effectiveLength => dst%Lambda_sl(:,en) - prm%w) - associate(BoltzmannRatio => prm%Q_s/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & - effectiveLength => dst%Lambda_sl(:,en) - prm%w) + significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) + StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) + vel = prm%h/(t_n + t_k) - vel = prm%h/(t_n + t_k) + dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal + else where significantPositiveTau + dot_gamma_pos = 0.0_pReal + end where significantPositiveTau - dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal - else where significantPositiveTau - dot_gamma_pos = 0.0_pReal - end where significantPositiveTau + if (present(ddot_gamma_dtau_pos)) then + significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) + dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls + dtk = -1.0_pReal * t_k / tau_pos - if (present(ddot_gamma_dtau_pos)) then - significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) - dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls - dtk = -1.0_pReal * t_k / tau_pos + dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal + ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal + else where significantPositiveTau2 + ddot_gamma_dtau_pos = 0.0_pReal + end where significantPositiveTau2 + end if - ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal - else where significantPositiveTau2 - ddot_gamma_dtau_pos = 0.0_pReal - end where significantPositiveTau2 - endif + significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) + StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) + vel = prm%h/(t_n + t_k) - vel = prm%h/(t_n + t_k) + dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal + else where significantNegativeTau + dot_gamma_neg = 0.0_pReal + end where significantNegativeTau - dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal - else where significantNegativeTau - dot_gamma_neg = 0.0_pReal - end where significantNegativeTau + if (present(ddot_gamma_dtau_neg)) then + significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) + dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & + * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls + dtk = -1.0_pReal * t_k / tau_neg - if (present(ddot_gamma_dtau_neg)) then - significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) - dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls - dtk = -1.0_pReal * t_k / tau_neg + dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal - dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal + ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal + else where significantNegativeTau2 + ddot_gamma_dtau_neg = 0.0_pReal + end where significantNegativeTau2 + end if - ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal - else where significantNegativeTau2 - ddot_gamma_dtau_neg = 0.0_pReal - end where significantNegativeTau2 - end if - - end associate + end associate end associate end subroutine kinetics From bff186051c47e1f62170cb234a5a4fa704c4c8b4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 30 Oct 2021 22:08:55 +0200 Subject: [PATCH 05/28] simplified --- src/phase_mechanical_plastic_dislotungsten.f90 | 16 +++++++--------- src/phase_mechanical_plastic_dislotwin.f90 | 3 ++- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index fce013bb2..6b2e65796 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -462,7 +462,7 @@ pure subroutine kinetics(Mp,T,ph,en, & real(pReal), dimension(param(ph)%sum_N_sl) :: & StressRatio, & StressRatio_p,StressRatio_pminus1, & - dvel, vel, & + dvel, & tau_pos,tau_neg, & t_n, t_k, dtk,dtn integer :: j @@ -487,12 +487,11 @@ pure subroutine kinetics(Mp,T,ph,en, & StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) + t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) & + / (prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - vel = prm%h/(t_n + t_k) - - dot_gamma_pos = dot_gamma_0 * sign(vel,tau_pos) * 0.5_pReal + dot_gamma_pos = dot_gamma_0 * sign(prm%h/(t_n + t_k),tau_pos) * 0.5_pReal else where significantPositiveTau dot_gamma_pos = 0.0_pReal end where significantPositiveTau @@ -516,12 +515,11 @@ pure subroutine kinetics(Mp,T,ph,en, & StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) - t_n = prm%b_sl/(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)*prm%omega*effectiveLength) + t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) & + / (prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) - vel = prm%h/(t_n + t_k) - - dot_gamma_neg = dot_gamma_0 * sign(vel,tau_neg) * 0.5_pReal + dot_gamma_neg = dot_gamma_0 * sign(prm%h/(t_n + t_k),tau_neg) * 0.5_pReal else where significantNegativeTau dot_gamma_neg = 0.0_pReal end where significantNegativeTau diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index ac179d775..de73cee04 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -890,7 +890,8 @@ pure subroutine kinetics_sl(Mp,T,ph,en, & stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p Q_kB_T = prm%Q_sl/(kB*T) - v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) + v_wait_inverse = exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q) & + / prm%v_0 v_run_inverse = prm%B/(tau_eff*prm%b_sl) dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) From e20b705f5416ee81935b8b2e70747f33076ea2c5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 30 Oct 2021 23:02:51 +0200 Subject: [PATCH 06/28] following dislotwin --- src/phase_mechanical_plastic_dislotungsten.f90 | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 6b2e65796..6a947cc70 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -463,7 +463,7 @@ pure subroutine kinetics(Mp,T,ph,en, & StressRatio, & StressRatio_p,StressRatio_pminus1, & dvel, & - tau_pos,tau_neg, & + tau_pos, tau_neg, tau_eff, & t_n, t_k, dtk,dtn integer :: j @@ -482,8 +482,10 @@ pure subroutine kinetics(Mp,T,ph,en, & dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & effectiveLength => dst%Lambda_sl(:,en) - prm%w) - significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) - StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls + tau_eff = abs(tau_pos)-dst%tau_pass(:,en) + + significantPositiveTau: where(tau_eff > tol_math_check) + StressRatio = tau_eff/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) @@ -510,8 +512,10 @@ pure subroutine kinetics(Mp,T,ph,en, & end where significantPositiveTau2 end if - significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) - StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls + tau_eff = abs(tau_neg)-dst%tau_pass(:,en) + + significantNegativeTau: where(tau_eff > tol_math_check) + StressRatio = tau_eff/prm%tau_Peierls StressRatio_p = StressRatio** prm%p StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) From f5fe0b9dca43ff022f9de1154814821a5821b32d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 30 Oct 2021 23:03:27 +0200 Subject: [PATCH 07/28] bugfix: change of behavior negative values for the resolved stress do not make sense. The paper does not take this into account (eq (14), Cereceda et. al 2016). According to my understanding, only the non-thermal contributions should be substracted, so abs(tau_pos)/abs(tau_neg) would not be sufficient. --- src/phase_mechanical_plastic_dislotungsten.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 6a947cc70..9f565df67 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -491,7 +491,7 @@ pure subroutine kinetics(Mp,T,ph,en, & t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) & / (prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff) ! corrected eq. (14) dot_gamma_pos = dot_gamma_0 * sign(prm%h/(t_n + t_k),tau_pos) * 0.5_pReal else where significantPositiveTau @@ -521,7 +521,7 @@ pure subroutine kinetics(Mp,T,ph,en, & t_n = prm%b_sl*exp(BoltzmannRatio*(1.0_pReal-StressRatio_p) ** prm%q) & / (prm%omega*effectiveLength) - t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_pos) + t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff) ! corrected eq. (14) dot_gamma_neg = dot_gamma_0 * sign(prm%h/(t_n + t_k),tau_neg) * 0.5_pReal else where significantNegativeTau From 25ca77c38e96a6e99f0c1fc295efc4d77e1dc4e2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 13:32:24 +0100 Subject: [PATCH 08/28] parameters from original paper (mostly) --- examples/config/Phase_Dislotungsten_W.yaml | 26 ---------------- .../mechanical/plastic/dislotungsten_W.yaml | 30 +++++++++++++++++++ 2 files changed, 30 insertions(+), 26 deletions(-) delete mode 100644 examples/config/Phase_Dislotungsten_W.yaml create mode 100644 examples/config/phase/mechanical/plastic/dislotungsten_W.yaml diff --git a/examples/config/Phase_Dislotungsten_W.yaml b/examples/config/Phase_Dislotungsten_W.yaml deleted file mode 100644 index bf8796cfa..000000000 --- a/examples/config/Phase_Dislotungsten_W.yaml +++ /dev/null @@ -1,26 +0,0 @@ -type: dislotungsten - -N_sl: [12] - -rho_mob_0: [1.0e+9] -rho_dip_0: [1.0] - -nu_a: [9.1e+11] -b_sl: [2.72e-10] -Delta_H_kp,0: [2.61154e-19] # 1.63 eV, Delta_H0 - -tau_Peierls: [2.03e+9] -p_sl: [0.86] -q_sl: [1.69] -h: [2.566e-10] -w: [2.992e-09] -B: [8.3e-5] -D_a: 1.0 # d_edge - -# climb (disabled) -D_0: 0.0 -Q_cl: 0.0 -V_cl: [0.0] - -h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] -a_nonSchmid: [0.938, 0.71, 4.43] diff --git a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml new file mode 100644 index 000000000..ea6ea5383 --- /dev/null +++ b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml @@ -0,0 +1,30 @@ +type: dislotungsten +references: + - D. Cereceda et al., + International Journal of Plasticity 78:242-265, 2016, + http://dx.doi.org/10.1016/j.ijplas.2015.09.002 +N_sl: [12] +b_sl: [2.72e-10] +rho_mob_0: [1.0e+9] # estimated from section 3.2 +rho_dip_0: [1.0] # not given +Q_s: [2.61154e-19] # 1.63 eV, Delta_H0 +B: [8.3e-5] +omega: [9.1e+11] # nu_0 +p_sl: [0.86] +q_sl: [1.69] +tau_Peierls: [2.03e+9] # there is also tau_c^* = 2.92GPa, not clear in the paper +h: [2.566e-10] +a_nonSchmid: [0.938, 0.71, 4.43] +h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] +w: [2.992e-09] # 11b +i_sl: [1] # c, eq. (25) +D: 1.0e+10 # d_g, eq. (25) +D_a: 1.0 # d_edge = D_a*b + +# climb (disabled) +D_0: 0.0 # disable climb +f_at: 1 +Q_cl: 1.0 +output: [Lambda_sl] + +v_0: [1] From fbc4865c30c2a306fdec8af0f11c8fcfdaba22ad Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 15:42:18 +0100 Subject: [PATCH 09/28] mathematically equivalent re-formulation allows to disable contribution from grain size by setting it to a large value --- .../config/phase/mechanical/plastic/dislotungsten_W.yaml | 3 +-- src/phase_mechanical_plastic_dislotungsten.f90 | 9 +++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml index ea6ea5383..6a0e4dcc2 100644 --- a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml +++ b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml @@ -18,10 +18,9 @@ a_nonSchmid: [0.938, 0.71, 4.43] h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] w: [2.992e-09] # 11b i_sl: [1] # c, eq. (25) -D: 1.0e+10 # d_g, eq. (25) +D: 1.0e+20 # d_g, eq. (25) D_a: 1.0 # d_edge = D_a*b -# climb (disabled) D_0: 0.0 # disable climb f_at: 1 Q_cl: 1.0 diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 9f565df67..21ea34611 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -196,7 +196,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl ! sanity checks - if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' + if ( prm%D_0 < 0.0_pReal) extmsg = trim(extmsg)//' D_0' if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' @@ -375,16 +375,17 @@ module subroutine dislotungsten_dependentState(ph,en) en real(pReal), dimension(param(ph)%sum_N_sl) :: & - dislocationSpacing + Lambda_sl_inv associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) - dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%tau_pass(:,en) = prm%mu*prm%b_sl & * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) - dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl) + Lambda_sl_inv = 1.0_pReal/prm%D & + + sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl + dst%Lambda_sl(:,en) = Lambda_sl_inv**(-1.0_pReal) end associate From a352f8deebcef3a8c4681066a55c8dc625dd9358 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 16:05:22 +0100 Subject: [PATCH 10/28] reasonable convergence --- .../phase/mechanical/plastic/dislotungsten_W.yaml | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml index 6a0e4dcc2..6efe77644 100644 --- a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml +++ b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml @@ -3,6 +3,10 @@ references: - D. Cereceda et al., International Journal of Plasticity 78:242-265, 2016, http://dx.doi.org/10.1016/j.ijplas.2015.09.002 + - A. Koester et al., + Acta Materialia 60:3894-3901, 2012 + http://dx.doi.org/10.1016/j.actamat.2012.03.053 +output: [Lambda_sl] N_sl: [12] b_sl: [2.72e-10] rho_mob_0: [1.0e+9] # estimated from section 3.2 @@ -12,18 +16,19 @@ B: [8.3e-5] omega: [9.1e+11] # nu_0 p_sl: [0.86] q_sl: [1.69] -tau_Peierls: [2.03e+9] # there is also tau_c^* = 2.92GPa, not clear in the paper +tau_Peierls: [2.03e+9] h: [2.566e-10] -a_nonSchmid: [0.938, 0.71, 4.43] h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] w: [2.992e-09] # 11b i_sl: [1] # c, eq. (25) D: 1.0e+20 # d_g, eq. (25) D_a: 1.0 # d_edge = D_a*b +# values in Cereceda et al. are high, using parameters from Koester et al. +a_nonSchmid: [0.61, 0.23, 0.55] # Table 1 + D_0: 0.0 # disable climb f_at: 1 Q_cl: 1.0 -output: [Lambda_sl] v_0: [1] From 6d117882c5965e0b5cc9583bd46bce5779b17a46 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 17:29:35 +0100 Subject: [PATCH 11/28] use precalculated value --- src/phase_mechanical_plastic_dislotungsten.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 21ea34611..d1b932619 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -502,7 +502,7 @@ pure subroutine kinetics(Mp,T,ph,en, & if (present(ddot_gamma_dtau_pos)) then significantPositiveTau2: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls + * StressRatio_pminus1 / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_pos dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal @@ -532,7 +532,7 @@ pure subroutine kinetics(Mp,T,ph,en, & if (present(ddot_gamma_dtau_neg)) then significantNegativeTau2: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check) dtn = -1.0_pReal * t_n * BoltzmannRatio * prm%p * prm%q * (1.0_pReal-StressRatio_p)**(prm%q - 1.0_pReal) & - * (StressRatio)**(prm%p - 1.0_pReal) / prm%tau_Peierls + * StressRatio_pminus1 / prm%tau_Peierls dtk = -1.0_pReal * t_k / tau_neg dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal From 737eff91793e8b794124878f5a22bfb42b0699e4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 17:31:01 +0100 Subject: [PATCH 12/28] bugfix: change of behavior v_0 was erroneously introduced use tau = (tau_pos+tau_neg)/2 --- PRIVATE | 2 +- .../phase/mechanical/plastic/dislotungsten_W.yaml | 2 -- src/phase_mechanical_plastic_dislotungsten.f90 | 12 ++++-------- 3 files changed, 5 insertions(+), 11 deletions(-) diff --git a/PRIVATE b/PRIVATE index fabe69749..00a536a78 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit fabe69749425e8a7aceb3b7c2758b40d97d8b809 +Subproject commit 00a536a78508cb273071517128a7edc7c387088b diff --git a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml index 6efe77644..8f474d2a3 100644 --- a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml +++ b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml @@ -30,5 +30,3 @@ a_nonSchmid: [0.61, 0.23, 0.55] # Table 1 D_0: 0.0 # disable climb f_at: 1 Q_cl: 1.0 - -v_0: [1] diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index d1b932619..8970fd53c 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -24,7 +24,6 @@ submodule(phase:plastic) dislotungsten tau_Peierls, & !< Peierls stress !* mobility law parameters Q_s, & !< activation energy for glide [J] - v_0, & !< dislocation velocity prefactor [m/s] p, & !< p-exponent in glide velocity q, & !< q-exponent in glide velocity B, & !< friction coefficient @@ -158,7 +157,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl)) rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl)) - prm%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl)) prm%b_sl = pl%get_as1dFloat('b_sl', requiredSize=size(N_sl)) prm%Q_s = pl%get_as1dFloat('Q_s', requiredSize=size(N_sl)) @@ -189,7 +187,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%w = math_expand(prm%w, N_sl) prm%omega = math_expand(prm%omega, N_sl) prm%tau_Peierls = math_expand(prm%tau_Peierls, N_sl) - prm%v_0 = math_expand(prm%v_0, N_sl) prm%B = math_expand(prm%B, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl) prm%f_at = math_expand(prm%f_at, N_sl) @@ -200,7 +197,6 @@ module function plastic_dislotungsten_init() result(myPlasticity) if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' - if (any(prm%v_0 < 0.0_pReal)) extmsg = trim(extmsg)//' v_0' if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' if (any(prm%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s' if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls' @@ -211,7 +207,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) else slipActive rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray allocate(prm%b_sl,prm%d_caron,prm%i_sl,prm%f_at,prm%tau_Peierls, & - prm%Q_s,prm%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & + prm%Q_s,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, & source = emptyRealArray) allocate(prm%forestProjection(0,0)) allocate(prm%h_sl_sl (0,0)) @@ -338,11 +334,11 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en) dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg) - where(dEq0(tau_pos)) ! ToDo: use avg of +/- + where(dEq0((tau_pos+tau_neg)*0.5_pReal)) dot_rho_dip_formation = 0.0_pReal dot_rho_dip_climb = 0.0_pReal else where - d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & ! ToDo: use avg of +/- + d_hat = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos+tau_neg)*0.5_pReal), & prm%d_caron, & ! lower limit dst%Lambda_sl(:,en)) ! upper limit dot_rho_dip_formation = merge(2.0_pReal*(d_hat-prm%d_caron)*stt%rho_mob(:,en)*dot%gamma_sl(:,en)/prm%b_sl, & @@ -480,7 +476,7 @@ pure subroutine kinetics(Mp,T,ph,en, & if (present(tau_neg_out)) tau_neg_out = tau_neg associate(BoltzmannRatio => prm%Q_s/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & + dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl, & effectiveLength => dst%Lambda_sl(:,en) - prm%w) tau_eff = abs(tau_pos)-dst%tau_pass(:,en) From 3f3224a9cb0734df1de18f04211f40d590dc43ec Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 18:59:57 +0100 Subject: [PATCH 13/28] found better source --- .../phase/mechanical/plastic/dislotungsten_W.yaml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml index 8f474d2a3..cb1d837a1 100644 --- a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml +++ b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml @@ -3,9 +3,9 @@ references: - D. Cereceda et al., International Journal of Plasticity 78:242-265, 2016, http://dx.doi.org/10.1016/j.ijplas.2015.09.002 - - A. Koester et al., - Acta Materialia 60:3894-3901, 2012 - http://dx.doi.org/10.1016/j.actamat.2012.03.053 + - R. Gröger et al., + Acta Materialia 56(19):5412-5425, 2008, + https://doi.org/10.1016/j.actamat.2008.07.037 output: [Lambda_sl] N_sl: [12] b_sl: [2.72e-10] @@ -24,8 +24,8 @@ i_sl: [1] # c, eq. (25) D: 1.0e+20 # d_g, eq. (25) D_a: 1.0 # d_edge = D_a*b -# values in Cereceda et al. are high, using parameters from Koester et al. -a_nonSchmid: [0.61, 0.23, 0.55] # Table 1 +# values in Cereceda et al. are high, using parameters from Gröger et al. +a_nonSchmid: [0.0, 0.56, 0.75] # Table 2 D_0: 0.0 # disable climb f_at: 1 From e0b6a28b4800a2438baa19aa9f8285df203c654b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 22:13:36 +0100 Subject: [PATCH 14/28] better matching name --- src/phase_mechanical_plastic_dislotungsten.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 8970fd53c..c759cdaad 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -476,7 +476,7 @@ pure subroutine kinetics(Mp,T,ph,en, & if (present(tau_neg_out)) tau_neg_out = tau_neg associate(BoltzmannRatio => prm%Q_s/(kB*T), & - dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl, & + b_rho_half => stt%rho_mob(:,en) * prm%b_sl * 0.5_pReal, & effectiveLength => dst%Lambda_sl(:,en) - prm%w) tau_eff = abs(tau_pos)-dst%tau_pass(:,en) @@ -490,7 +490,7 @@ pure subroutine kinetics(Mp,T,ph,en, & / (prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff) ! corrected eq. (14) - dot_gamma_pos = dot_gamma_0 * sign(prm%h/(t_n + t_k),tau_pos) * 0.5_pReal + dot_gamma_pos = b_rho_half * sign(prm%h/(t_n + t_k),tau_pos) else where significantPositiveTau dot_gamma_pos = 0.0_pReal end where significantPositiveTau @@ -503,7 +503,7 @@ pure subroutine kinetics(Mp,T,ph,en, & dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal - ddot_gamma_dtau_pos = dot_gamma_0 * dvel* 0.5_pReal + ddot_gamma_dtau_pos = b_rho_half * dvel else where significantPositiveTau2 ddot_gamma_dtau_pos = 0.0_pReal end where significantPositiveTau2 @@ -520,7 +520,7 @@ pure subroutine kinetics(Mp,T,ph,en, & / (prm%omega*effectiveLength) t_k = effectiveLength * prm%B /(2.0_pReal*prm%b_sl*tau_eff) ! corrected eq. (14) - dot_gamma_neg = dot_gamma_0 * sign(prm%h/(t_n + t_k),tau_neg) * 0.5_pReal + dot_gamma_neg = b_rho_half * sign(prm%h/(t_n + t_k),tau_neg) else where significantNegativeTau dot_gamma_neg = 0.0_pReal end where significantNegativeTau @@ -533,7 +533,7 @@ pure subroutine kinetics(Mp,T,ph,en, & dvel = -1.0_pReal * prm%h * (dtk + dtn) / (t_n + t_k)**2.0_pReal - ddot_gamma_dtau_neg = dot_gamma_0 * dvel * 0.5_pReal + ddot_gamma_dtau_neg = b_rho_half * dvel else where significantNegativeTau2 ddot_gamma_dtau_neg = 0.0_pReal end where significantNegativeTau2 From 509835bf0bed72a2e3842e368caeb80b996d2e95 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 22:37:54 +0100 Subject: [PATCH 15/28] type hints for tensor adjustments to Rotation/Orientation needed to enable type checking with mypy --- python/damask/_orientation.py | 2 +- python/damask/_rotation.py | 9 +++++++-- python/damask/tensor.py | 12 ++++++------ 3 files changed, 14 insertions(+), 9 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 3d4d259ff..bb14fd38b 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -125,7 +125,7 @@ class Orientation(Rotation,Crystal): """Create deep copy.""" dup = copy.deepcopy(self) if rotation is not None: - dup.quaternion = Orientation(rotation,family='cubic').quaternion + dup.quaternion = Rotation(rotation).quaternion return dup copy = __copy__ diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 74a3f7419..ac921d70a 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,3 +1,5 @@ +import copy + import numpy as np from . import tensor @@ -85,9 +87,12 @@ class Rotation: + str(self.quaternion) - def __copy__(self,**kwargs): + def __copy__(self,rotation=None): """Create deep copy.""" - return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion) + dup = copy.deepcopy(self) + if rotation is not None: + dup.quaternion = Rotation(rotation).quaternion + return dup copy = __copy__ diff --git a/python/damask/tensor.py b/python/damask/tensor.py index cf5d94020..a735b355e 100644 --- a/python/damask/tensor.py +++ b/python/damask/tensor.py @@ -8,7 +8,7 @@ All routines operate on numpy.ndarrays of shape (...,3,3). import numpy as _np -def deviatoric(T): +def deviatoric(T: _np.ndarray) -> _np.ndarray: """ Calculate deviatoric part of a tensor. @@ -26,7 +26,7 @@ def deviatoric(T): return T - spherical(T,tensor=True) -def eigenvalues(T_sym): +def eigenvalues(T_sym: _np.ndarray) -> _np.ndarray: """ Eigenvalues, i.e. principal components, of a symmetric tensor. @@ -45,7 +45,7 @@ def eigenvalues(T_sym): return _np.linalg.eigvalsh(symmetric(T_sym)) -def eigenvectors(T_sym,RHS=False): +def eigenvectors(T_sym: _np.ndarray, RHS: bool = False) -> _np.ndarray: """ Eigenvectors of a symmetric tensor. @@ -70,7 +70,7 @@ def eigenvectors(T_sym,RHS=False): return v -def spherical(T,tensor=True): +def spherical(T: _np.ndarray, tensor: bool = True) -> _np.ndarray: """ Calculate spherical part of a tensor. @@ -92,7 +92,7 @@ def spherical(T,tensor=True): return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph -def symmetric(T): +def symmetric(T: _np.ndarray) -> _np.ndarray: """ Symmetrize tensor. @@ -110,7 +110,7 @@ def symmetric(T): return (T+transpose(T))*0.5 -def transpose(T): +def transpose(T: _np.ndarray) -> _np.ndarray: """ Transpose tensor. From 2d25dfcdf2811d1671e66d57bfc94fb5da2f8f9b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 22:43:06 +0100 Subject: [PATCH 16/28] automated type checking --- .gitlab-ci.yml | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 8680b2873..1fd5a6555 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -81,7 +81,7 @@ checkout: - release ################################################################################################### -processing: +pytest: stage: python script: - cd $DAMASKROOT/python @@ -91,6 +91,16 @@ processing: - master - release +mypy: + stage: python + script: + - cd $DAMASKROOT/python + - mypy damask/mechanics.py + except: + - master + - release + + ################################################################################################### compile_grid_Intel: From 0bc267c76b6101e9886101b522343eff6445df41 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 31 Oct 2021 22:50:41 +0100 Subject: [PATCH 17/28] automated type checking for mechanics --- .gitlab-ci.yml | 3 +-- python/damask/mechanics.py | 42 ++++++++++++++++++++------------------ python/mypy.ini | 16 +++++++++++++++ 3 files changed, 39 insertions(+), 22 deletions(-) create mode 100644 python/mypy.ini diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 1fd5a6555..83e9f8934 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -95,13 +95,12 @@ mypy: stage: python script: - cd $DAMASKROOT/python - - mypy damask/mechanics.py + - mypy damask/tensor.py damask/mechanics.py except: - master - release - ################################################################################################### compile_grid_Intel: stage: compile diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 0e160523b..1a03f390b 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -5,13 +5,15 @@ All routines operate on numpy.ndarrays of shape (...,3,3). """ -from . import tensor as _tensor -from . import _rotation +from typing import Sequence import numpy as _np +from . import tensor as _tensor +from . import _rotation -def deformation_Cauchy_Green_left(F): + +def deformation_Cauchy_Green_left(F: _np.ndarray) -> _np.ndarray: """ Calculate left Cauchy-Green deformation tensor (Finger deformation tensor). @@ -29,7 +31,7 @@ def deformation_Cauchy_Green_left(F): return _np.matmul(F,_tensor.transpose(F)) -def deformation_Cauchy_Green_right(F): +def deformation_Cauchy_Green_right(F: _np.ndarray) -> _np.ndarray: """ Calculate right Cauchy-Green deformation tensor. @@ -47,7 +49,7 @@ def deformation_Cauchy_Green_right(F): return _np.matmul(_tensor.transpose(F),F) -def equivalent_strain_Mises(epsilon): +def equivalent_strain_Mises(epsilon: _np.ndarray) -> _np.ndarray: """ Calculate the Mises equivalent of a strain tensor. @@ -65,7 +67,7 @@ def equivalent_strain_Mises(epsilon): return _equivalent_Mises(epsilon,2.0/3.0) -def equivalent_stress_Mises(sigma): +def equivalent_stress_Mises(sigma: _np.ndarray) -> _np.ndarray: """ Calculate the Mises equivalent of a stress tensor. @@ -83,7 +85,7 @@ def equivalent_stress_Mises(sigma): return _equivalent_Mises(sigma,3.0/2.0) -def maximum_shear(T_sym): +def maximum_shear(T_sym: _np.ndarray) -> _np.ndarray: """ Calculate the maximum shear component of a symmetric tensor. @@ -102,7 +104,7 @@ def maximum_shear(T_sym): return (w[...,0] - w[...,2])*0.5 -def rotation(T): +def rotation(T: _np.ndarray) -> _rotation.Rotation: """ Calculate the rotational part of a tensor. @@ -120,7 +122,7 @@ def rotation(T): return _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0]) -def strain(F,t,m): +def strain(F: _np.ndarray, t: str, m: float) -> _np.ndarray: """ Calculate strain tensor (Seth–Hill family). @@ -160,7 +162,7 @@ def strain(F,t,m): return eps -def stress_Cauchy(P,F): +def stress_Cauchy(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Calculate the Cauchy stress (true stress). @@ -182,7 +184,7 @@ def stress_Cauchy(P,F): return _tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F)) -def stress_second_Piola_Kirchhoff(P,F): +def stress_second_Piola_Kirchhoff(P: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Calculate the second Piola-Kirchhoff stress. @@ -205,7 +207,7 @@ def stress_second_Piola_Kirchhoff(P,F): return _tensor.symmetric(_np.einsum('...ij,...jk',_np.linalg.inv(F),P)) -def stretch_left(T): +def stretch_left(T: _np.ndarray) -> _np.ndarray: """ Calculate left stretch of a tensor. @@ -223,7 +225,7 @@ def stretch_left(T): return _polar_decomposition(T,'V')[0] -def stretch_right(T): +def stretch_right(T: _np.ndarray) -> _np.ndarray: """ Calculate right stretch of a tensor. @@ -241,7 +243,7 @@ def stretch_right(T): return _polar_decomposition(T,'U')[0] -def _polar_decomposition(T,requested): +def _polar_decomposition(T: _np.ndarray, requested: Sequence[str]) -> tuple: """ Perform singular value decomposition. @@ -257,21 +259,21 @@ def _polar_decomposition(T,requested): u, _, vh = _np.linalg.svd(T) R = _np.einsum('...ij,...jk',u,vh) - output = () + output = [] if 'R' in requested: - output+=(R,) + output+=[R] if 'V' in requested: - output+=(_np.einsum('...ij,...kj',T,R),) + output+=[_np.einsum('...ij,...kj',T,R)] if 'U' in requested: - output+=(_np.einsum('...ji,...jk',R,T),) + output+=[_np.einsum('...ji,...jk',R,T)] if len(output) == 0: raise ValueError('output needs to be out of V, R, U') - return output + return tuple(output) -def _equivalent_Mises(T_sym,s): +def _equivalent_Mises(T_sym: _np.ndarray, s: float) -> _np.ndarray: """ Base equation for Mises equivalent of a stress or strain tensor. diff --git a/python/mypy.ini b/python/mypy.ini new file mode 100644 index 000000000..e6900781c --- /dev/null +++ b/python/mypy.ini @@ -0,0 +1,16 @@ +[mypy-scipy.*] +ignore_missing_imports = True +[mypy-h5py.*] +ignore_missing_imports = True +[mypy-vtk.*] +ignore_missing_imports = True +[mypy-PIL.*] +ignore_missing_imports = True +[mypy-matplotlib.*] +ignore_missing_imports = True +[mypy-yaml.*] +ignore_missing_imports = True +[mypy-pandas.*] +ignore_missing_imports = True +[mypy-wx.*] +ignore_missing_imports = True From 34e04fa45e98c522280c5832758d907c07d8445e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 1 Nov 2021 10:00:13 +0100 Subject: [PATCH 18/28] tests should work after updating test server typehints for numpy and pyyaml are now available --- python/mypy.ini | 2 -- 1 file changed, 2 deletions(-) diff --git a/python/mypy.ini b/python/mypy.ini index e6900781c..01001daa6 100644 --- a/python/mypy.ini +++ b/python/mypy.ini @@ -8,8 +8,6 @@ ignore_missing_imports = True ignore_missing_imports = True [mypy-matplotlib.*] ignore_missing_imports = True -[mypy-yaml.*] -ignore_missing_imports = True [mypy-pandas.*] ignore_missing_imports = True [mypy-wx.*] From 2e3de727ccf396d021116f2398d20126c2bbccfb Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 2 Nov 2021 07:27:08 +0100 Subject: [PATCH 19/28] hint at issues with the current parametrization --- .../phase/mechanical/plastic/dislotungsten_W.yaml | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml index cb1d837a1..55814c3f8 100644 --- a/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml +++ b/examples/config/phase/mechanical/plastic/dislotungsten_W.yaml @@ -20,13 +20,16 @@ tau_Peierls: [2.03e+9] h: [2.566e-10] h_sl-sl: [0.009, 0.72, 0.009, 0.05, 0.05, 0.06, 0.09] w: [2.992e-09] # 11b -i_sl: [1] # c, eq. (25) -D: 1.0e+20 # d_g, eq. (25) -D_a: 1.0 # d_edge = D_a*b # values in Cereceda et al. are high, using parameters from Gröger et al. a_nonSchmid: [0.0, 0.56, 0.75] # Table 2 -D_0: 0.0 # disable climb +# (almost) no annhilation, adjustment needed for simulations beyond the yield point +i_sl: [1] # c, eq. (25) +D: 1.0e+20 # d_g, eq. (25) +D_a: 1.0 # d_edge = D_a*b + +# disable climb (not discussed in Cereceda et al.) +D_0: 0.0 f_at: 1 Q_cl: 1.0 From 493a0969eb15837177e999add4d366377da6793e Mon Sep 17 00:00:00 2001 From: Test User Date: Tue, 2 Nov 2021 16:42:50 +0100 Subject: [PATCH 20/28] [skip ci] updated version information after successful test of v3.0.0-alpha5-45-g1a558db83 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 98903601b..8c8e85679 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-31-gddb25ad0e +v3.0.0-alpha5-45-g1a558db83 From 32aff9d96619c68b85652efdb69157b2bd099851 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 2 Nov 2021 13:01:32 -0400 Subject: [PATCH 21/28] added typehints to seeds.py --- .gitlab-ci.yml | 2 +- python/damask/seeds.py | 16 ++++++++++------ 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 83e9f8934..e064cb0af 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -95,7 +95,7 @@ mypy: stage: python script: - cd $DAMASKROOT/python - - mypy damask/tensor.py damask/mechanics.py + - mypy damask/tensor.py damask/mechanics.py damask/seeds.py except: - master - release diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 26aa3084b..2203c9495 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,13 +1,16 @@ """Functionality for generation of seed points for Voronoi or Laguerre tessellation.""" from scipy import spatial as _spatial +from typing import Sequence + import numpy as _np from . import util as _util from . import grid_filters as _grid_filters +from . import _grid -def from_random(size,N_seeds,cells=None,rng_seed=None): +def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_seed=None) -> _np.ndarray: """ Place seeds randomly in space. @@ -41,7 +44,8 @@ def from_random(size,N_seeds,cells=None,rng_seed=None): return coords -def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=None): +def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distance: float, + periodic: bool = True, rng_seed=None) -> _np.ndarray: """ Place seeds according to a Poisson disc distribution. @@ -75,18 +79,17 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed= i = 0 progress = _util._ProgressBar(N_seeds+1,'',50) while s < N_seeds: + i += 1 candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3)) tree = _spatial.cKDTree(coords[:s],boxsize=size) if periodic else \ _spatial.cKDTree(coords[:s]) distances = tree.query(candidates)[0] best = distances.argmax() if distances[best] > distance: # require minimum separation + i = 0 coords[s] = candidates[best] # maximum separation to existing point cloud s += 1 progress.update(s) - i = 0 - else: - i += 1 if i == 100: raise ValueError('Seeding not possible') @@ -94,7 +97,8 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed= return coords -def from_grid(grid,selection=None,invert=False,average=False,periodic=True): +def from_grid(grid: _grid.Grid, selection: Sequence[int] = None, + invert: bool = False, average: bool = False, periodic: bool = True) -> tuple[_np.ndarray, _np.ndarray]: """ Create seeds from grid description. From 735952bd32013b182b9c599c77e8d2b37ebd1404 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 2 Nov 2021 13:40:09 -0400 Subject: [PATCH 22/28] fixed typo --- python/damask/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 001e46276..231fa8b30 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -14,8 +14,8 @@ from . import tensor # noqa from . import mechanics # noqa from . import solver # noqa from . import grid_filters # noqa -#Modules that contain only one class (of the same name), are prefixed by a '_'. -#For example, '_colormap' containsa class called 'Colormap' which is imported as 'damask.Colormap'. +# Modules that contain only one class (of the same name), are prefixed by a '_'. +# For example, '_colormap' contains a class called 'Colormap' which is imported as 'damask.Colormap'. from ._rotation import Rotation # noqa from ._crystal import Crystal # noqa from ._orientation import Orientation # noqa From 8636c5dad4be71fee5805c6394634ee7b146bbeb Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 2 Nov 2021 13:41:05 -0400 Subject: [PATCH 23/28] removed grid typehint <-- no clue how to break circular import --- python/damask/seeds.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 2203c9495..85b28850b 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -7,7 +7,6 @@ import numpy as _np from . import util as _util from . import grid_filters as _grid_filters -from . import _grid def from_random(size: _np.ndarray, N_seeds: int, cells: _np.ndarray = None, rng_seed=None) -> _np.ndarray: @@ -97,7 +96,7 @@ def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distan return coords -def from_grid(grid: _grid.Grid, selection: Sequence[int] = None, +def from_grid(grid, selection: Sequence[int] = None, invert: bool = False, average: bool = False, periodic: bool = True) -> tuple[_np.ndarray, _np.ndarray]: """ Create seeds from grid description. @@ -105,15 +104,15 @@ def from_grid(grid: _grid.Grid, selection: Sequence[int] = None, Parameters ---------- grid : damask.Grid - Grid, from which the material IDs are used as seeds. + Grid from which the material IDs are used as seeds. selection : iterable of integers, optional Material IDs to consider. invert : boolean, false - Do not consider the material IDs given in selection. Defaults to False. + Consider all material IDs except those in selection. Defaults to False. average : boolean, optional Seed corresponds to center of gravity of material ID cloud. periodic : boolean, optional - Center of gravity with periodic boundaries. + Center of gravity accounts for periodic boundaries. Returns ------- From ccfe276ae1283a2806b0b1dd9ad08119d7371b7e Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 2 Nov 2021 14:28:54 -0400 Subject: [PATCH 24/28] fixed type hinting for seeds.py and grid_filters.py --- python/damask/grid_filters.py | 56 ++++++++++++++++++++--------------- python/damask/seeds.py | 4 +-- 2 files changed, 34 insertions(+), 26 deletions(-) diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index 816c727cd..cb6388f41 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -12,9 +12,11 @@ the following operations are required for tensorial data: """ from scipy import spatial as _spatial +from typing import Sequence, Tuple, Union + import numpy as _np -def _ks(size,cells,first_order=False): +def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order: bool = False) -> _np.ndarray: """ Get wave numbers operator. @@ -41,7 +43,7 @@ def _ks(size,cells,first_order=False): return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1) -def curl(size,f): +def curl(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray: u""" Calculate curl of a vector or tensor field in Fourier space. @@ -72,7 +74,7 @@ def curl(size,f): return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3]) -def divergence(size,f): +def divergence(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray: u""" Calculate divergence of a vector or tensor field in Fourier space. @@ -99,7 +101,7 @@ def divergence(size,f): return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3]) -def gradient(size,f): +def gradient(size: _np.ndarray, f: _np.ndarray) -> _np.ndarray: u""" Calculate gradient of a scalar or vector fieldin Fourier space. @@ -126,7 +128,9 @@ def gradient(size,f): return _np.fft.irfftn(grad_,axes=(0,1,2),s=f.shape[:3]) -def coordinates0_point(cells,size,origin=_np.zeros(3)): +def coordinates0_point(cells: Union[ _np.ndarray,Sequence[int]], + size: _np.ndarray, + origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Cell center positions (undeformed). @@ -145,8 +149,8 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)): Undeformed cell center coordinates. """ - start = origin + size/cells*.5 - end = origin + size - size/cells*.5 + start = origin + size/_np.array(cells)*.5 + end = origin + size - size/_np.array(cells)*.5 return _np.stack(_np.meshgrid(_np.linspace(start[0],end[0],cells[0]), _np.linspace(start[1],end[1],cells[1]), @@ -154,7 +158,7 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)): axis = -1) -def displacement_fluct_point(size,F): +def displacement_fluct_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Cell center displacement field from fluctuation part of the deformation gradient field. @@ -186,7 +190,7 @@ def displacement_fluct_point(size,F): return _np.fft.irfftn(displacement,axes=(0,1,2),s=F.shape[:3]) -def displacement_avg_point(size,F): +def displacement_avg_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Cell center displacement field from average part of the deformation gradient field. @@ -207,7 +211,7 @@ def displacement_avg_point(size,F): return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size)) -def displacement_point(size,F): +def displacement_point(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Cell center displacement field from deformation gradient field. @@ -227,7 +231,7 @@ def displacement_point(size,F): return displacement_avg_point(size,F) + displacement_fluct_point(size,F) -def coordinates_point(size,F,origin=_np.zeros(3)): +def coordinates_point(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Cell center positions. @@ -249,7 +253,8 @@ def coordinates_point(size,F,origin=_np.zeros(3)): return coordinates0_point(F.shape[:3],size,origin) + displacement_point(size,F) -def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True): +def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray, + ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]: """ Return grid 'DNA', i.e. cells, size, and origin from 1D array of point positions. @@ -292,13 +297,15 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True): raise ValueError('Regular cell spacing violated.') if ordered and not _np.allclose(coordinates0.reshape(tuple(cells)+(3,),order='F'), - coordinates0_point(cells,size,origin),atol=atol): + coordinates0_point(list(cells),size,origin),atol=atol): raise ValueError('Input data is not ordered (x fast, z slow).') return (cells,size,origin) -def coordinates0_node(cells,size,origin=_np.zeros(3)): +def coordinates0_node(cells: Union[_np.ndarray,Sequence[int]], + size: _np.ndarray, + origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Nodal positions (undeformed). @@ -323,7 +330,7 @@ def coordinates0_node(cells,size,origin=_np.zeros(3)): axis = -1) -def displacement_fluct_node(size,F): +def displacement_fluct_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Nodal displacement field from fluctuation part of the deformation gradient field. @@ -343,7 +350,7 @@ def displacement_fluct_node(size,F): return point_to_node(displacement_fluct_point(size,F)) -def displacement_avg_node(size,F): +def displacement_avg_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Nodal displacement field from average part of the deformation gradient field. @@ -364,7 +371,7 @@ def displacement_avg_node(size,F): return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size)) -def displacement_node(size,F): +def displacement_node(size: _np.ndarray, F: _np.ndarray) -> _np.ndarray: """ Nodal displacement field from deformation gradient field. @@ -384,7 +391,7 @@ def displacement_node(size,F): return displacement_avg_node(size,F) + displacement_fluct_node(size,F) -def coordinates_node(size,F,origin=_np.zeros(3)): +def coordinates_node(size: _np.ndarray, F: _np.ndarray, origin: _np.ndarray = _np.zeros(3)) -> _np.ndarray: """ Nodal positions. @@ -406,7 +413,8 @@ def coordinates_node(size,F,origin=_np.zeros(3)): return coordinates0_node(F.shape[:3],size,origin) + displacement_node(size,F) -def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True): +def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray, + ordered: bool = True) -> Tuple[_np.ndarray,_np.ndarray,_np.ndarray]: """ Return grid 'DNA', i.e. cells, size, and origin from 1D array of nodal positions. @@ -441,13 +449,13 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True): raise ValueError('Regular cell spacing violated.') if ordered and not _np.allclose(coordinates0.reshape(tuple(cells+1)+(3,),order='F'), - coordinates0_node(cells,size,origin),atol=atol): + coordinates0_node(list(cells),size,origin),atol=atol): raise ValueError('Input data is not ordered (x fast, z slow).') return (cells,size,origin) -def point_to_node(cell_data): +def point_to_node(cell_data: _np.ndarray) -> _np.ndarray: """ Interpolate periodic point data to nodal data. @@ -469,7 +477,7 @@ def point_to_node(cell_data): return _np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap') -def node_to_point(node_data): +def node_to_point(node_data: _np.ndarray) -> _np.ndarray: """ Interpolate periodic nodal data to point data. @@ -491,7 +499,7 @@ def node_to_point(node_data): return c[1:,1:,1:] -def coordinates0_valid(coordinates0): +def coordinates0_valid(coordinates0: _np.ndarray) -> bool: """ Check whether coordinates form a regular grid. @@ -513,7 +521,7 @@ def coordinates0_valid(coordinates0): return False -def regrid(size,F,cells): +def regrid(size: _np.ndarray, F: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]]) -> _np.ndarray: """ Return mapping from coordinates in deformed configuration to a regular grid. diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 85b28850b..1d7dfc58e 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,7 +1,7 @@ """Functionality for generation of seed points for Voronoi or Laguerre tessellation.""" from scipy import spatial as _spatial -from typing import Sequence +from typing import Sequence,Tuple import numpy as _np @@ -97,7 +97,7 @@ def from_Poisson_disc(size: _np.ndarray, N_seeds: int, N_candidates: int, distan def from_grid(grid, selection: Sequence[int] = None, - invert: bool = False, average: bool = False, periodic: bool = True) -> tuple[_np.ndarray, _np.ndarray]: + invert: bool = False, average: bool = False, periodic: bool = True) -> Tuple[_np.ndarray, _np.ndarray]: """ Create seeds from grid description. From 59a6dc365277d9eb3b997c5cd43f15bf18034ba4 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 2 Nov 2021 14:40:09 -0400 Subject: [PATCH 25/28] run mypy on all python/damask --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index e064cb0af..14783263d 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -95,7 +95,7 @@ mypy: stage: python script: - cd $DAMASKROOT/python - - mypy damask/tensor.py damask/mechanics.py damask/seeds.py + - mypy damask/*.py except: - master - release From f0c587d4aa468446201cdf850735a98728547332 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 3 Nov 2021 07:53:38 +0100 Subject: [PATCH 26/28] polishing oder of imports is build-in, 3rd party, internal --- .gitlab-ci.yml | 2 +- python/damask/grid_filters.py | 4 +++- python/damask/seeds.py | 2 +- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 14783263d..165cdc68f 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -95,7 +95,7 @@ mypy: stage: python script: - cd $DAMASKROOT/python - - mypy damask/*.py + - mypy -m damask except: - master - release diff --git a/python/damask/grid_filters.py b/python/damask/grid_filters.py index cb6388f41..42b5a16c4 100644 --- a/python/damask/grid_filters.py +++ b/python/damask/grid_filters.py @@ -11,11 +11,13 @@ the following operations are required for tensorial data: - D1 = D3.reshape(cells+(-1,)).reshape(-1,9,order='F') """ -from scipy import spatial as _spatial + from typing import Sequence, Tuple, Union +from scipy import spatial as _spatial import numpy as _np + def _ks(size: _np.ndarray, cells: Union[_np.ndarray,Sequence[int]], first_order: bool = False) -> _np.ndarray: """ Get wave numbers operator. diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 1d7dfc58e..4d5a8c624 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,8 +1,8 @@ """Functionality for generation of seed points for Voronoi or Laguerre tessellation.""" -from scipy import spatial as _spatial from typing import Sequence,Tuple +from scipy import spatial as _spatial import numpy as _np from . import util as _util From ba66d7b8164693d31efb94d36d87e1e5678710bf Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 4 Nov 2021 21:19:55 +0100 Subject: [PATCH 27/28] [skip ci] updated version information after successful test of v3.0.0-alpha5-52-g5b0ad1eb8 --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 8c8e85679..9ddbae920 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-45-g1a558db83 +v3.0.0-alpha5-52-g5b0ad1eb8 From 5142813a9879117cd11652bab75c10b8bbd55ba8 Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 5 Nov 2021 10:23:36 +0100 Subject: [PATCH 28/28] [skip ci] updated version information after successful test of v3.0.0-alpha5-64-g8e08af31e --- python/damask/VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/VERSION b/python/damask/VERSION index 9ddbae920..bd451b119 100644 --- a/python/damask/VERSION +++ b/python/damask/VERSION @@ -1 +1 @@ -v3.0.0-alpha5-52-g5b0ad1eb8 +v3.0.0-alpha5-64-g8e08af31e