223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge remote-tracking branch 'origin/development' into do-concurrent

This commit is contained in:
Martin Diehl 2021-07-26 20:20:52 +02:00
parent d19ab4c4f6 2be24351e5
commit c7ed0e7934
22 changed files with 565 additions and 606 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit d92b030e5777e718c77edc2e1e93abfa0981b024 Subproject commit 72c58103860e127d37ccf3a06827331de29406ca

42
examples/config/Phase_DisloUCLA_Tungsten.config
View File

@ -1,42 +0,0 @@
[Tungsten]
elasticity hooke
plasticity disloucla
(output) edge_density
(output) dipole_density
(output) shear_rate_slip
(output) accumulated_shear_slip
(output) resolved_stress_slip
(output) threshold_stress_slip
grainsize 2.7e-5 # Average grain size [m] 2.0e-5
SolidSolutionStrength 0.0 # Strength due to elements in solid solution
### Dislocation glide parameters ###
#per family
Nslip 12
slipburgers 2.72e-10 # Burgers vector of slip system [m]
rhoedge0 1.0e12 # Initial edge dislocation density [m/m**3]
rhoedgedip0 1.0 # Initial edged dipole dislocation density [m/m**3]
Qedge 2.61154e-19 # Activation energy for dislocation glide [J], 1.63 eV
v0 1 # Initial glide velocity [m/s]
p_slip 0.86 # p-exponent in glide velocity
q_slip 1.69 # q-exponent in glide velocity
tau_peierls 2.03e9 # peierls stress [Pa]
#mobility law
kink_height 2.567e-10 # kink height sqrt(6)/3*lattice_parameter [m]
omega 9.1e11 # attemp frequency (from kMC paper) [s^(-1)]
kink_width 29.95e-10 # kink pair width ~ 11 b (kMC paper) [m]
dislolength 78e-10 # dislocation length (ideally lambda) [m] initial value 11b
friction_coeff 8.3e-5 # friction coeff. B [Pa*s]
#hardening
dipoleformationfactor 0 # to have hardening due to dipole formation off
CLambdaSlip 10.0 # Adj. parameter controlling dislocation mean free path
D0 4.0e-5 # Vacancy diffusion prefactor [m**2/s]
Qsd 4.5e-19 # Activation energy for climb [J]
Catomicvolume 1.0 # Adj. parameter controlling the atomic volume [in b]
Cedgedipmindistance 1.0 # Adj. parameter controlling the minimum dipole distance [in b]
interaction_slipslip 0.009 0.72 0.009 0.05 0.05 0.06 0.09
nonschmid_coefficients 0.938 0.71 4.43 0.0 0.0 0.0

26
examples/config/Phase_Dislotungsten_W.yaml Normal file
View File

@ -0,0 +1,26 @@
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]

2
examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml
View File

@ -6,7 +6,7 @@ b_sl: [2.56e-10]
rho_mob_0: [1.0e+12] rho_mob_0: [1.0e+12]
rho_dip_0: [1.0] rho_dip_0: [1.0]
v_0: [1.0e+4] v_0: [1.0e+4]
Q_s: [3.7e-19] Q_sl: [3.7e-19]
p_sl: [1.0] p_sl: [1.0]
q_sl: [1.0] q_sl: [1.0]
tau_0: [1.5e+8] tau_0: [1.5e+8]

3
examples/config/phase/mechanical/plastic/dislotwin_IF-steel.yaml
View File

@ -11,7 +11,7 @@ b_sl: [2.49e-10, 2.49e-10]
rho_mob_0: [2.81e12, 2.8e+12] rho_mob_0: [2.81e12, 2.8e+12]
rho_dip_0: [1.0, 1.0] # not given rho_dip_0: [1.0, 1.0] # not given
v_0: [1.4e+3, 1.4e+3] v_0: [1.4e+3, 1.4e+3]
Q_s: [1.57e-19, 1.57e-19] # Delta_F Q_sl: [1.57e-19, 1.57e-19] # Delta_F
tau_0: [454.e+6, 454.e+6] tau_0: [454.e+6, 454.e+6]
p_sl: [0.325, 0.325] p_sl: [0.325, 0.325]
q_sl: [1.55, 1.55] q_sl: [1.55, 1.55]
@ -19,6 +19,5 @@ i_sl: [23.3, 23.3]
D_a: 7.4 # C_anni D_a: 7.4 # C_anni
B: [0.001, 0.001] B: [0.001, 0.001]
h_sl-sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073] h_sl-sl: [0.1, 0.72, 0.1, 0.053, 0.053, 0.073, 0.137, 0.72, 0.72, 0.053, 0.053, 0.053, 0.053, 0.073, 0.073, 0.073, 0.073, 0.073, 0.073, 0.137, 0.073, 0.073, 0.137, 0.073]
D_0: 4.0e-05
Q_cl: 5.4e-19 # no recovery! Q_cl: 5.4e-19 # no recovery!
D: 40.e-6 # estimated D: 40.e-6 # estimated

33
examples/config/Phase_Phenopowerlaw_Magnesium.yaml → examples/config/phase/mechanical/plastic/phenopowerlaw_Mg.yaml
View File

@ -1,20 +1,29 @@
N_sl: [3, 3, 0, 6, 0, 6]
N_tw: [6, 0, 0, 6]
h_0_tw-tw: 50.0e+6
h_0_sl-sl: 500.0e+6
h_0_tw-sl: 150.0e+6
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
output: [xi_sl, xi_tw]
type: phenopowerlaw type: phenopowerlaw
xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6] references:
- F. Wang et al.,
Acta Materialia 80:77-93, 2014,
https://doi.org/10.1016/j.actamat.2014.07.048
output: [xi_sl, xi_tw]
N_sl: [3, 3, 0, 6, 0, 6] # basal, 1. prism, -, 1. pyr<a>, -, 2. pyr<c+a>
N_tw: [6, 0, 6] # tension, -, compression
xi_0_sl: [10.e+6, 55.e+6, 0., 60.e+6, 0., 60.e+6]
xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6] xi_inf_sl: [40.e+6, 135.e+6, 0., 150.e+6, 0., 150.e+6]
xi_0_tw: [40.e+6, 0., 0., 60.e+6] xi_0_tw: [40.e+6, 0., 60.e+6]
a_sl: 2.25 a_sl: 2.25
dot_gamma_0_sl: 0.001 dot_gamma_0_sl: 0.001
dot_gamma_0_tw: 0.001 dot_gamma_0_tw: 0.001
n_sl: 20 n_sl: 20
n_tw: 20 n_tw: 20
f_sat_sl-tw: 10.0 f_sat_sl-tw: 10.0
h_0_sl-sl: 500.0e+6
h_0_tw-tw: 50.0e+6
h_0_tw-sl: 150.0e+6
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_tw-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
h_sl-tw: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

5
examples/config/phase/mechanical/plastic/phenopowerlaw_Ti.yaml
View File

@ -7,14 +7,17 @@ references:
Acta Materialia 132:598-610, 2017, Acta Materialia 132:598-610, 2017,
https://doi.org/10.1016/j.actamat.2017.05.015 https://doi.org/10.1016/j.actamat.2017.05.015
output: [gamma_sl] output: [gamma_sl]
N_sl: [3, 3, 0, 0, 12]
N_sl: [3, 3, 0, 0, 12] # basal, 1. prism, -, -, 2. pyr<c+a>
n_sl: 20 n_sl: 20
a_sl: 2.0 a_sl: 2.0
dot_gamma_0_sl: 0.001 dot_gamma_0_sl: 0.001
h_0_sl-sl: 200.e+6 h_0_sl-sl: 200.e+6
# C. Zambaldi et al.: # C. Zambaldi et al.:
xi_0_sl: [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6] xi_0_sl: [349.e+6, 150.e+6, 0.0, 0.0, 1107.e+6]
xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6] xi_inf_sl: [568.e+6, 150.e+7, 0.0, 0.0, 3420.e+6]
# L. Wang et al. : # L. Wang et al. :
# xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6] # xi_0_sl: [127.e+6, 96.e+6, 0.0, 0.0, 240.e+6]
h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] h_sl-sl: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]

4
examples/mesh/tensionZ.load
View File

@ -1,11 +1,11 @@
#initial elastic step #initial elastic step
$Loadcase 1 time 0.0005 incs 1 frequency 5 $Loadcase 1 t 0.0005 N 1 f_out 5
Face 1 X 0.01 Face 1 X 0.01
Face 2 X 0.0 Face 2 X 0.0
Face 2 Y 0.0 Face 2 Y 0.0
Face 2 Z 0.0 Face 2 Z 0.0
$EndLoadcase $EndLoadcase
$Loadcase 2 time 10.0 incs 200 frequency 5 $Loadcase 2 t 10.0 N 200 f_out 5
Face 1 X 0.01 Face 1 X 0.01
Face 2 X 0.0 Face 2 X 0.0
Face 2 Y 0.0 Face 2 Y 0.0

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha4-137-gb69b85754 v3.0.0-alpha4-182-gac6d31b1f

9
src/YAML_types.f90
View File

@ -849,7 +849,7 @@ function tNode_get_byKey_as1dFloat(self,k,defaultVal,requiredSize) result(nodeAs
if (self%contains(k)) then if (self%contains(k)) then
node => self%get(k) node => self%get(k)
select type(self) select type(node)
class is(tList) class is(tList)
list => node%asList() list => node%asList()
nodeAs1dFloat = list%as1dFloat() nodeAs1dFloat = list%as1dFloat()
@ -872,11 +872,12 @@ end function tNode_get_byKey_as1dFloat
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Access by key and convert to float array (2D) !> @brief Access by key and convert to float array (2D)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat) function tNode_get_byKey_as2dFloat(self,k,defaultVal,requiredShape) result(nodeAs2dFloat)
class(tNode), intent(in), target :: self class(tNode), intent(in), target :: self
character(len=*), intent(in) :: k character(len=*), intent(in) :: k
real(pReal), intent(in), dimension(:,:), optional :: defaultVal real(pReal), intent(in), dimension(:,:), optional :: defaultVal
integer, intent(in), dimension(2), optional :: requiredShape
real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat real(pReal), dimension(:,:), allocatable :: nodeAs2dFloat
@ -898,6 +899,10 @@ function tNode_get_byKey_as2dFloat(self,k,defaultVal) result(nodeAs2dFloat)
call IO_error(143,ext_msg=k) call IO_error(143,ext_msg=k)
endif endif
if (present(requiredShape)) then
if (any(requiredShape /= shape(nodeAs2dFloat))) call IO_error(146,ext_msg=k)
endif
end function tNode_get_byKey_as2dFloat end function tNode_get_byKey_as2dFloat

13
src/math.f90
View File

@ -545,19 +545,6 @@ pure function math_symmetric33(m)
end function math_symmetric33 end function math_symmetric33
!--------------------------------------------------------------------------------------------------
!> @brief symmetrize a 6x6 matrix
!--------------------------------------------------------------------------------------------------
pure function math_symmetric66(m)
real(pReal), dimension(6,6) :: math_symmetric66
real(pReal), dimension(6,6), intent(in) :: m
math_symmetric66 = 0.5_pReal * (m + transpose(m))
end function math_symmetric66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief skew part of a 3x3 matrix !> @brief skew part of a 3x3 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

178
src/mesh/DAMASK_mesh.f90
View File

@ -22,6 +22,15 @@ program DAMASK_mesh
implicit none implicit none
type :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1 !< frequency of result writes
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file ! variables related to information from load case and geom file
integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
@ -68,7 +77,7 @@ program DAMASK_mesh
type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases type(tLoadCase), allocatable, dimension(:) :: loadCases !< array of all load cases
type(tSolutionState), allocatable, dimension(:) :: solres type(tSolutionState), allocatable, dimension(:) :: solres
PetscInt :: faceSet, currentFaceSet, field, dimPlex PetscInt :: faceSet, currentFaceSet, dimPlex
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: & external :: &
@ -91,8 +100,7 @@ program DAMASK_mesh
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D) call DMGetDimension(geomMesh,dimPlex,ierr) !< dimension of mesh (2D or 3D)
CHKERRA(ierr) CHKERRA(ierr)
nActiveFields = 1 allocate(solres(1))
allocate(solres(nActiveFields))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! reading basic information from load case file and allocate data structure containing load cases ! reading basic information from load case file and allocate data structure containing load cases
@ -103,8 +111,8 @@ program DAMASK_mesh
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_stringValue(line,chunkPos,i))
case('$loadcase') case('$Loadcase')
N_def = N_def + 1 N_def = N_def + 1
end select end select
enddo ! count all identifiers to allocate memory and do sanity check enddo ! count all identifiers to allocate memory and do sanity check
@ -114,32 +122,26 @@ program DAMASK_mesh
allocate(loadCases(N_def)) allocate(loadCases(N_def))
do i = 1, size(loadCases) do i = 1, size(loadCases)
allocate(loadCases(i)%fieldBC(nActiveFields)) allocate(loadCases(i)%fieldBC(1))
field = 1 loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID
loadCases(i)%fieldBC(field)%ID = FIELD_MECH_ID
enddo enddo
do i = 1, size(loadCases) do i = 1, size(loadCases)
do field = 1, nActiveFields loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements
select case (loadCases(i)%fieldBC(field)%ID) allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents))
case(FIELD_MECH_ID) do component = 1, loadCases(i)%fieldBC(1)%nComponents
loadCases(i)%fieldBC(field)%nComponents = dimPlex !< X, Y (, Z) displacements select case (component)
allocate(loadCases(i)%fieldBC(field)%componentBC(loadCases(i)%fieldBC(field)%nComponents)) case (1)
do component = 1, loadCases(i)%fieldBC(field)%nComponents loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID
select case (component) case (2)
case (1) loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_X_ID case (3)
case (2) loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3)
loadCases(i)%fieldBC(field)%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select
enddo
end select end select
do component = 1, loadCases(i)%fieldBC(field)%nComponents enddo
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal) do component = 1, loadCases(i)%fieldBC(1)%nComponents
allocate(loadCases(i)%fieldBC(field)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pReal)
enddo allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
enddo enddo
enddo enddo
@ -151,52 +153,45 @@ program DAMASK_mesh
chunkPos = IO_stringPos(line) chunkPos = IO_stringPos(line)
do i = 1, chunkPos(1) do i = 1, chunkPos(1)
select case (IO_lc(IO_stringValue(line,chunkPos,i))) select case (IO_stringValue(line,chunkPos,i))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loadcase information ! loadcase information
case('$loadcase') case('$Loadcase')
currentLoadCase = IO_intValue(line,chunkPos,i+1) currentLoadCase = IO_intValue(line,chunkPos,i+1)
case('face') case('Face')
currentFace = IO_intValue(line,chunkPos,i+1) currentFace = IO_intValue(line,chunkPos,i+1)
currentFaceSet = -1 currentFaceSet = -1
do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
enddo enddo
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC') if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
case('t','time','delta') ! increment time case('t')
loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1) loadCases(currentLoadCase)%time = IO_floatValue(line,chunkPos,i+1)
case('n','incs','increments','steps') ! number of increments case('N')
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1) loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) case('f_out')
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
loadCases(currentLoadCase)%logscale = 1
case('freq','frequency','outputfreq') ! frequency of result writings
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1) loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('guessreset','dropguessing') case('estimate_rate')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! boundary condition information ! boundary condition information
case('x','y','z') case('X','Y','Z')
select case(IO_lc(IO_stringValue(line,chunkPos,i))) select case(IO_stringValue(line,chunkPos,i))
case('x') case('X')
ID = COMPONENT_MECH_X_ID ID = COMPONENT_MECH_X_ID
case('y') case('Y')
ID = COMPONENT_MECH_Y_ID ID = COMPONENT_MECH_Y_ID
case('z') case('Z')
ID = COMPONENT_MECH_Z_ID ID = COMPONENT_MECH_Z_ID
end select end select
do field = 1, nActiveFields do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents
if (loadCases(currentLoadCase)%fieldBC(field)%ID == FIELD_MECH_ID) then if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then
do component = 1, loadcases(currentLoadCase)%fieldBC(field)%nComponents loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = &
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%ID == ID) then .true.
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask (currentFaceSet) = & loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = &
.true. IO_floatValue(line,chunkPos,i+1)
loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Value(currentFaceSet) = &
IO_floatValue(line,chunkPos,i+1)
endif
enddo
endif endif
enddo enddo
end select end select
@ -212,21 +207,16 @@ program DAMASK_mesh
print'(a,i0)', ' load case: ', currentLoadCase print'(a,i0)', ' load case: ', currentLoadCase
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(currentLoadCase)%followFormerTrajectory) &
print'(a)', ' drop guessing along trajectory' print'(a)', ' drop guessing along trajectory'
do field = 1, nActiveFields print'(a)', ' Field '//trim(FIELD_MECH_label)
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
print'(a)', ' Field '//trim(FIELD_MECH_label)
end select do faceSet = 1, mesh_Nboundaries
do faceSet = 1, mesh_Nboundaries do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents
do component = 1, loadCases(currentLoadCase)%fieldBC(field)%nComponents if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) &
if (loadCases(currentLoadCase)%fieldBC(field)%componentBC(component)%Mask(faceSet)) & print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), &
print'(a,i2,a,i2,a,f12.7)', ' Face ', mesh_boundaries(faceSet), & ' Component ', component, &
' Component ', component, & ' Value ', loadCases(currentLoadCase)%fieldBC(1)% &
' Value ', loadCases(currentLoadCase)%fieldBC(field)% & componentBC(component)%Value(faceSet)
componentBC(component)%Value(faceSet) enddo
enddo
enddo
enddo enddo
print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time print'(a,f12.6)', ' time: ', loadCases(currentLoadCase)%time
if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count
@ -240,12 +230,7 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers ! doing initialization depending on active solvers
call FEM_Utilities_init call FEM_Utilities_init
do field = 1, nActiveFields call FEM_mechanical_init(loadCases(1)%fieldBC(1))
select case (loadCases(1)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mechanical_init(loadCases(1)%fieldBC(field))
end select
enddo
if (worldrank == 0) then if (worldrank == 0) then
open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE') open(newunit=statUnit,file=trim(getSolverJobName())//'.sta',form='FORMATTED',status='REPLACE')
@ -266,32 +251,14 @@ program DAMASK_mesh
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end timeIncOld = timeinc ! last timeinc that brought former inc to an end
if (loadCases(currentLoadCase)%logscale == 0) then ! linear scale timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal)
else
if (currentLoadCase == 1) then ! 1st load case of logarithmic scale
if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd
else ! not-1st inc of 1st load case of logarithmic scale
timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal))
endif
else ! not-1st load case of logarithmic scale
timeinc = time0 * &
( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal))&
-(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc-1 ,pReal)/&
real(loadCases(currentLoadCase)%incs ,pReal)))
endif
endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time
time = time + timeinc ! forward target time time = time + timeinc ! forward target time
stepFraction = stepFraction + 1 ! count step stepFraction = stepFraction + 1 ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
@ -306,33 +273,16 @@ program DAMASK_mesh
'-',stepFraction, '/', subStepFactor**cutBackLevel '-',stepFraction, '/', subStepFactor**cutBackLevel
flush(IO_STDOUT) flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------- call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
! forward fields
do field = 1, nActiveFields
select case (loadCases(currentLoadCase)%fieldBC(field)%ID)
case(FIELD_MECH_ID)
call FEM_mechanical_forward (&
guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve fields ! solve fields
stagIter = 0 stagIter = 0
stagIterate = .true. stagIterate = .true.
do while (stagIterate) do while (stagIterate)
do field = 1, nActiveFields solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1))
select case (loadCases(currentLoadCase)%fieldBC(field)%ID) if(.not. solres(1)%converged) exit
case(FIELD_MECH_ID)
solres(field) = FEM_mechanical_solution (&
incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(field))
end select
if(.not. solres(field)%converged) exit ! no solution found
enddo
stagIter = stagIter + 1 stagIter = stagIter + 1
stagIterate = stagIter < stagItMax & stagIterate = stagIter < stagItMax &
.and. all(solres(:)%converged) & .and. all(solres(:)%converged) &

340
src/mesh/FEM_quadrature.f90
View File

@ -9,7 +9,7 @@ module FEM_quadrature
private private
integer, parameter :: & integer, parameter :: &
maxOrder = 5 !< current max interpolation set at cubic (intended to be arbitrary) maxOrder = 5 !< maximum integration order
real(pReal), dimension(2,3), parameter :: & real(pReal), dimension(2,3), parameter :: &
triangle = reshape([-1.0_pReal, -1.0_pReal, & triangle = reshape([-1.0_pReal, -1.0_pReal, &
1.0_pReal, -1.0_pReal, & 1.0_pReal, -1.0_pReal, &
@ -20,8 +20,12 @@ module FEM_quadrature
-1.0_pReal, 1.0_pReal, -1.0_pReal, & -1.0_pReal, 1.0_pReal, -1.0_pReal, &
-1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4]) -1.0_pReal, -1.0_pReal, 1.0_pReal], shape=[3,4])
type :: group_float !< variable length datatype
real(pReal), dimension(:), allocatable :: p
end type group_float
integer, dimension(2:3,maxOrder), public, protected :: & integer, dimension(2:3,maxOrder), public, protected :: &
FEM_nQuadrature !< number of quadrature points for a given spatial dimension(2-3) and interpolation order(1-maxOrder) FEM_nQuadrature !< number of quadrature points for spatial dimension(2-3) and interpolation order (1-maxOrder)
type(group_float), dimension(2:3,maxOrder), public, protected :: & type(group_float), dimension(2:3,maxOrder), public, protected :: &
FEM_quadrature_weights, & !< quadrature weights for each quadrature rule FEM_quadrature_weights, & !< quadrature weights for each quadrature rule
FEM_quadrature_points !< quadrature point coordinates (in simplical system) for each quadrature rule FEM_quadrature_points !< quadrature point coordinates (in simplical system) for each quadrature rule
@ -35,145 +39,146 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief initializes FEM interpolation data !> @brief initializes FEM interpolation data
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_quadrature_init subroutine FEM_quadrature_init()
print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6) print'(/,a)', ' <<<+- FEM_quadrature init -+>>>'; flush(6)
print*, 'L. Zhang et al., Journal of Computational Mathematics 27(1):89-96, 2009'
print*, 'https://www.jstor.org/stable/43693493'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D linear ! 2D linear
FEM_nQuadrature(2,1) = 1 FEM_nQuadrature(2,1) = 1
allocate(FEM_quadrature_weights(2,1)%p(1)) allocate(FEM_quadrature_weights(2,1)%p(FEM_nQuadrature(2,1)))
FEM_quadrature_weights(2,1)%p(1) = 1.0_pReal FEM_quadrature_weights(2,1)%p(1) = 1._pReal
allocate(FEM_quadrature_points (2,1)%p(2)) FEM_quadrature_points (2,1)%p = permutationStar3([1._pReal/3._pReal])
FEM_quadrature_points (2,1)%p(1:2) = permutationStar3([1.0_pReal/3.0_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D quadratic ! 2D quadratic
FEM_nQuadrature(2,2) = 3 FEM_nQuadrature(2,2) = 3
allocate(FEM_quadrature_weights(2,2)%p(3)) allocate(FEM_quadrature_weights(2,2)%p(FEM_nQuadrature(2,2)))
FEM_quadrature_weights(2,2)%p(1:3) = 1.0_pReal/3.0_pReal FEM_quadrature_weights(2,2)%p(1:3) = 1._pReal/3._pReal
allocate(FEM_quadrature_points (2,2)%p(6)) FEM_quadrature_points (2,2)%p = permutationStar21([1._pReal/6._pReal])
FEM_quadrature_points (2,2)%p(1:6) = permutationStar21([1.0_pReal/6.0_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D cubic ! 2D cubic
FEM_nQuadrature(2,3) = 6 FEM_nQuadrature(2,3) = 6
allocate(FEM_quadrature_weights(2,3)%p(6)) allocate(FEM_quadrature_weights(2,3)%p(FEM_nQuadrature(2,3)))
FEM_quadrature_weights(2,3)%p(1:3) = 0.22338158967801146570_pReal FEM_quadrature_weights(2,3)%p(1:3) = 2.2338158967801147e-1_pReal
FEM_quadrature_weights(2,3)%p(4:6) = 0.10995174365532186764_pReal FEM_quadrature_weights(2,3)%p(4:6) = 1.0995174365532187e-1_pReal
allocate(FEM_quadrature_points (2,3)%p(12)) FEM_quadrature_points (2,3)%p = [ &
FEM_quadrature_points (2,3)%p(1:6) = permutationStar21([0.44594849091596488632_pReal]) permutationStar21([4.4594849091596489e-1_pReal]), &
FEM_quadrature_points (2,3)%p(7:12)= permutationStar21([0.091576213509770743460_pReal]) permutationStar21([9.157621350977074e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D quartic ! 2D quartic
FEM_nQuadrature(2,4) = 12 FEM_nQuadrature(2,4) = 12
allocate(FEM_quadrature_weights(2,4)%p(12)) allocate(FEM_quadrature_weights(2,4)%p(FEM_nQuadrature(2,4)))
FEM_quadrature_weights(2,4)%p(1:3) = 0.11678627572638_pReal FEM_quadrature_weights(2,4)%p(1:3) = 1.1678627572637937e-1_pReal
FEM_quadrature_weights(2,4)%p(4:6) = 0.05084490637021_pReal FEM_quadrature_weights(2,4)%p(4:6) = 5.0844906370206817e-2_pReal
FEM_quadrature_weights(2,4)%p(7:12) = 0.08285107561837_pReal FEM_quadrature_weights(2,4)%p(7:12) = 8.285107561837358e-2_pReal
allocate(FEM_quadrature_points (2,4)%p(24)) FEM_quadrature_points (2,4)%p = [ &
FEM_quadrature_points (2,4)%p(1:6) = permutationStar21([0.24928674517091_pReal]) permutationStar21([2.4928674517091042e-1_pReal]), &
FEM_quadrature_points (2,4)%p(7:12) = permutationStar21([0.06308901449150_pReal]) permutationStar21([6.308901449150223e-2_pReal]), &
FEM_quadrature_points (2,4)%p(13:24)= permutationStar111([0.31035245103378_pReal, 0.63650249912140_pReal]) permutationStar111([3.1035245103378440e-1_pReal, 5.3145049844816947e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 2D quintic ! 2D quintic
FEM_nQuadrature(2,5) = 16 FEM_nQuadrature(2,5) = 16
allocate(FEM_quadrature_weights(2,5)%p(16)) allocate(FEM_quadrature_weights(2,5)%p(FEM_nQuadrature(2,5)))
FEM_quadrature_weights(2,5)%p(1 ) = 0.14431560767779_pReal FEM_quadrature_weights(2,5)%p(1:1) = 1.4431560767778717e-1_pReal
FEM_quadrature_weights(2,5)%p(2:4) = 0.09509163426728_pReal FEM_quadrature_weights(2,5)%p(2:4) = 9.509163426728463e-2_pReal
FEM_quadrature_weights(2,5)%p(5:7) = 0.10321737053472_pReal FEM_quadrature_weights(2,5)%p(5:7) = 1.0321737053471825e-1_pReal
FEM_quadrature_weights(2,5)%p(8:10) = 0.03245849762320_pReal FEM_quadrature_weights(2,5)%p(8:10) = 3.2458497623198080e-2_pReal
FEM_quadrature_weights(2,5)%p(11:16)= 0.02723031417443_pReal FEM_quadrature_weights(2,5)%p(11:16) = 2.7230314174434994e-2_pReal
allocate(FEM_quadrature_points (2,5)%p(32)) FEM_quadrature_points (2,5)%p = [ &
FEM_quadrature_points (2,5)%p(1:2) = permutationStar3([0.33333333333333_pReal]) permutationStar3([1._pReal/3._pReal]), &
FEM_quadrature_points (2,5)%p(3:8) = permutationStar21([0.45929258829272_pReal]) permutationStar21([4.5929258829272316e-1_pReal]), &
FEM_quadrature_points (2,5)%p(9:14) = permutationStar21([0.17056930775176_pReal]) permutationStar21([1.705693077517602e-1_pReal]), &
FEM_quadrature_points (2,5)%p(15:20)= permutationStar21([0.05054722831703_pReal]) permutationStar21([5.0547228317030975e-2_pReal]), &
FEM_quadrature_points (2,5)%p(21:32)= permutationStar111([0.26311282963464_pReal, 0.72849239295540_pReal]) permutationStar111([2.631128296346381e-1_pReal, 8.3947774099576053e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D linear ! 3D linear
FEM_nQuadrature(3,1) = 1 FEM_nQuadrature(3,1) = 1
allocate(FEM_quadrature_weights(3,1)%p(1)) allocate(FEM_quadrature_weights(3,1)%p(FEM_nQuadrature(3,1)))
FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal FEM_quadrature_weights(3,1)%p(1) = 1.0_pReal
allocate(FEM_quadrature_points (3,1)%p(3)) FEM_quadrature_points (3,1)%p = permutationStar4([0.25_pReal])
FEM_quadrature_points (3,1)%p(1:3)= permutationStar4([0.25_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D quadratic ! 3D quadratic
FEM_nQuadrature(3,2) = 4 FEM_nQuadrature(3,2) = 4
allocate(FEM_quadrature_weights(3,2)%p(4)) allocate(FEM_quadrature_weights(3,2)%p(FEM_nQuadrature(3,2)))
FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal FEM_quadrature_weights(3,2)%p(1:4) = 0.25_pReal
allocate(FEM_quadrature_points (3,2)%p(12)) FEM_quadrature_points (3,2)%p = permutationStar31([1.3819660112501052e-1_pReal])
FEM_quadrature_points (3,2)%p(1:12)= permutationStar31([0.13819660112501051518_pReal])
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D cubic ! 3D cubic
FEM_nQuadrature(3,3) = 14 FEM_nQuadrature(3,3) = 14
allocate(FEM_quadrature_weights(3,3)%p(14)) allocate(FEM_quadrature_weights(3,3)%p(FEM_nQuadrature(3,3)))
FEM_quadrature_weights(3,3)%p(5:8) = 0.11268792571801585080_pReal FEM_quadrature_weights(3,3)%p(1:4) = 7.3493043116361949e-2_pReal
FEM_quadrature_weights(3,3)%p(1:4) = 0.073493043116361949544_pReal FEM_quadrature_weights(3,3)%p(5:8) = 1.1268792571801585e-1_pReal
FEM_quadrature_weights(3,3)%p(9:14) = 0.042546020777081466438_pReal FEM_quadrature_weights(3,3)%p(9:14) = 4.2546020777081467e-2_pReal
allocate(FEM_quadrature_points (3,3)%p(42)) FEM_quadrature_points (3,3)%p = [ &
FEM_quadrature_points (3,3)%p(1:12) = permutationStar31([0.092735250310891226402_pReal]) permutationStar31([9.273525031089123e-2_pReal]), &
FEM_quadrature_points (3,3)%p(13:24)= permutationStar31([0.31088591926330060980_pReal]) permutationStar31([3.108859192633006e-1_pReal]), &
FEM_quadrature_points (3,3)%p(25:42)= permutationStar22([0.045503704125649649492_pReal]) permutationStar22([4.5503704125649649e-2_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D quartic ! 3D quartic (lower precision/unknown source)
FEM_nQuadrature(3,4) = 35 FEM_nQuadrature(3,4) = 35
allocate(FEM_quadrature_weights(3,4)%p(35)) allocate(FEM_quadrature_weights(3,4)%p(FEM_nQuadrature(3,4)))
FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal FEM_quadrature_weights(3,4)%p(1:4) = 0.0021900463965388_pReal
FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal FEM_quadrature_weights(3,4)%p(5:16) = 0.0143395670177665_pReal
FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal FEM_quadrature_weights(3,4)%p(17:22) = 0.0250305395686746_pReal
FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal FEM_quadrature_weights(3,4)%p(23:34) = 0.0479839333057554_pReal
FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal FEM_quadrature_weights(3,4)%p(35) = 0.0931745731195340_pReal
allocate(FEM_quadrature_points (3,4)%p(105)) FEM_quadrature_points (3,4)%p = [ &
FEM_quadrature_points (3,4)%p(1:12) = permutationStar31([0.0267367755543735_pReal]) permutationStar31([0.0267367755543735_pReal]), &
FEM_quadrature_points (3,4)%p(13:48) = permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]) permutationStar211([0.0391022406356488_pReal, 0.7477598884818090_pReal]), &
FEM_quadrature_points (3,4)%p(49:66) = permutationStar22([0.4547545999844830_pReal]) permutationStar22([0.4547545999844830_pReal]), &
FEM_quadrature_points (3,4)%p(67:102) = permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]) permutationStar211([0.2232010379623150_pReal, 0.0504792790607720_pReal]), &
FEM_quadrature_points (3,4)%p(103:105)= permutationStar4([0.25_pReal]) permutationStar4([0.25_pReal]) ]
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! 3D quintic ! 3D quintic (lower precision/unknown source)
FEM_nQuadrature(3,5) = 56 FEM_nQuadrature(3,5) = 56
allocate(FEM_quadrature_weights(3,5)%p(56)) allocate(FEM_quadrature_weights(3,5)%p(FEM_nQuadrature(3,5)))
FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal FEM_quadrature_weights(3,5)%p(1:4) = 0.0010373112336140_pReal
FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal FEM_quadrature_weights(3,5)%p(5:16) = 0.0096016645399480_pReal
FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal FEM_quadrature_weights(3,5)%p(17:28) = 0.0164493976798232_pReal
FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal FEM_quadrature_weights(3,5)%p(29:40) = 0.0153747766513310_pReal
FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal FEM_quadrature_weights(3,5)%p(41:52) = 0.0293520118375230_pReal
FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal FEM_quadrature_weights(3,5)%p(53:56) = 0.0366291366405108_pReal
allocate(FEM_quadrature_points (3,5)%p(168)) FEM_quadrature_points (3,5)%p = [ &
FEM_quadrature_points (3,5)%p(1:12) = permutationStar31([0.0149520651530592_pReal]) permutationStar31([0.0149520651530592_pReal]), &
FEM_quadrature_points (3,5)%p(13:48) = permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]) permutationStar211([0.0340960211962615_pReal, 0.1518319491659370_pReal]), &
FEM_quadrature_points (3,5)%p(49:84) = permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]) permutationStar211([0.0462051504150017_pReal, 0.3549340560639790_pReal]), &
FEM_quadrature_points (3,5)%p(85:120) = permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]) permutationStar211([0.2281904610687610_pReal, 0.0055147549744775_pReal]), &
FEM_quadrature_points (3,5)%p(121:156)= permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]) permutationStar211([0.3523052600879940_pReal, 0.0992057202494530_pReal]), &
FEM_quadrature_points (3,5)%p(157:168)= permutationStar31([0.1344783347929940_pReal]) permutationStar31([0.1344783347929940_pReal]) ]
call selfTest
end subroutine FEM_quadrature_init end subroutine FEM_quadrature_init
@ -186,11 +191,9 @@ pure function permutationStar3(point) result(qPt)
real(pReal), dimension(2) :: qPt real(pReal), dimension(2) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,1) :: temp
temp(:,1) = [point(1), point(1), point(1)] qPt = pack(matmul(triangle,reshape([ &
point(1), point(1), point(1)],[3,1])),.true.)
qPt = reshape(matmul(triangle, temp),[2])
end function permutationStar3 end function permutationStar3
@ -203,13 +206,11 @@ pure function permutationStar21(point) result(qPt)
real(pReal), dimension(6) :: qPt real(pReal), dimension(6) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(3,3) :: temp
temp(:,1) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1)] qPt = pack(matmul(triangle,reshape([ &
temp(:,2) = [point(1), 1.0_pReal - 2.0_pReal*point(1), point(1)] point(1), point(1), 1.0_pReal - 2.0_pReal*point(1), &
temp(:,3) = [1.0_pReal - 2.0_pReal*point(1), point(1), point(1)] point(1), 1.0_pReal - 2.0_pReal*point(1), point(1), &
1.0_pReal - 2.0_pReal*point(1), point(1), point(1)],[3,3])),.true.)
qPt = reshape(matmul(triangle, temp),[6])
end function permutationStar21 end function permutationStar21
@ -222,16 +223,14 @@ pure function permutationStar111(point) result(qPt)
real(pReal), dimension(12) :: qPt real(pReal), dimension(12) :: qPt
real(pReal), dimension(2), intent(in) :: point real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(3,6) :: temp
temp(:,1) = [point(1), point(2), 1.0_pReal - point(1) - point(2)] qPt = pack(matmul(triangle,reshape([ &
temp(:,2) = [point(1), 1.0_pReal - point(1) - point(2), point(2)] point(1), point(2), 1.0_pReal - point(1) - point(2), &
temp(:,3) = [point(2), point(1), 1.0_pReal - point(1) - point(2)] point(1), 1.0_pReal - point(1) - point(2), point(2), &
temp(:,4) = [point(2), 1.0_pReal - point(1) - point(2), point(1)] point(2), point(1), 1.0_pReal - point(1) - point(2), &
temp(:,5) = [1.0_pReal - point(1) - point(2), point(2), point(1)] point(2), 1.0_pReal - point(1) - point(2), point(1), &
temp(:,6) = [1.0_pReal - point(1) - point(2), point(1), point(2)] 1.0_pReal - point(1) - point(2), point(2), point(1), &
1.0_pReal - point(1) - point(2), point(1), point(2)],[3,6])),.true.)
qPt = reshape(matmul(triangle, temp),[12])
end function permutationStar111 end function permutationStar111
@ -244,11 +243,9 @@ pure function permutationStar4(point) result(qPt)
real(pReal), dimension(3) :: qPt real(pReal), dimension(3) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,1) :: temp
temp(:,1) = [point(1), point(1), point(1), point(1)] qPt = pack(matmul(tetrahedron,reshape([ &
point(1), point(1), point(1), point(1)],[4,1])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[3])
end function permutationStar4 end function permutationStar4
@ -261,14 +258,12 @@ pure function permutationStar31(point) result(qPt)
real(pReal), dimension(12) :: qPt real(pReal), dimension(12) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,4) :: temp
temp(:,1) = [point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2) = [point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1)] point(1), point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), &
temp(:,3) = [point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1)] point(1), point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), &
temp(:,4) = [1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)] point(1), 1.0_pReal - 3.0_pReal*point(1), point(1), point(1), &
1.0_pReal - 3.0_pReal*point(1), point(1), point(1), point(1)],[4,4])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[12])
end function permutationStar31 end function permutationStar31
@ -276,21 +271,19 @@ end function permutationStar31
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief star 22 permutation of input !> @brief star 22 permutation of input
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function permutationStar22(point) result(qPt) function permutationStar22(point) result(qPt)
real(pReal), dimension(18) :: qPt real(pReal), dimension(18) :: qPt
real(pReal), dimension(1), intent(in) :: point real(pReal), dimension(1), intent(in) :: point
real(pReal), dimension(4,6) :: temp
temp(:,1) = [point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2) = [point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1)] point(1), point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), &
temp(:,3) = [0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1)] point(1), 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), &
temp(:,4) = [0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1)] 0.5_pReal - point(1), point(1), point(1), 0.5_pReal - point(1), &
temp(:,5) = [0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1)] 0.5_pReal - point(1), point(1), 0.5_pReal - point(1), point(1), &
temp(:,6) = [point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)] 0.5_pReal - point(1), 0.5_pReal - point(1), point(1), point(1), &
point(1), 0.5_pReal - point(1), 0.5_pReal - point(1), point(1)],[4,6])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[18])
end function permutationStar22 end function permutationStar22
@ -303,22 +296,20 @@ pure function permutationStar211(point) result(qPt)
real(pReal), dimension(36) :: qPt real(pReal), dimension(36) :: qPt
real(pReal), dimension(2), intent(in) :: point real(pReal), dimension(2), intent(in) :: point
real(pReal), dimension(4,12) :: temp
temp(:,1 ) = [point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2 ) = [point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2)] point(1), point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
temp(:,3 ) = [point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), &
temp(:,4 ) = [point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] point(1), point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
temp(:,5 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2)] point(1), point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
temp(:,6 ) = [point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1)] point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), &
temp(:,7 ) = [point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2)] point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), &
temp(:,8 ) = [point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1)] point(2), point(1), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), &
temp(:,9 ) = [point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1)] point(2), point(1), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), &
temp(:,10) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2)] point(2), 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), &
temp(:,11) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1)] 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(1), point(2), &
temp(:,12) = [1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)] 1.0_pReal - 2.0_pReal*point(1) - point(2), point(1), point(2), point(1), &
1.0_pReal - 2.0_pReal*point(1) - point(2), point(2), point(1), point(1)],[4,12])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[36])
end function permutationStar211 end function permutationStar211
@ -331,35 +322,60 @@ pure function permutationStar1111(point) result(qPt)
real(pReal), dimension(72) :: qPt real(pReal), dimension(72) :: qPt
real(pReal), dimension(3), intent(in) :: point real(pReal), dimension(3), intent(in) :: point
real(pReal), dimension(4,24) :: temp
temp(:,1 ) = [point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3)] qPt = pack(matmul(tetrahedron,reshape([ &
temp(:,2 ) = [point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3)] point(1), point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,3 ) = [point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3)] point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
temp(:,4 ) = [point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2)] point(1), point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,5 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3)] point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
temp(:,6 ) = [point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2)] point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), &
temp(:,7 ) = [point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3)] point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), &
temp(:,8 ) = [point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3)] point(2), point(1), point(3), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,9 ) = [point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3)] point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(3), &
temp(:,10) = [point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1)] point(2), point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,11) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3)] point(2), point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
temp(:,12) = [point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1)] point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), &
temp(:,13) = [point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3)] point(2), 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), &
temp(:,14) = [point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2)] point(3), point(1), point(2), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,15) = [point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3)] point(3), point(1), 1.0_pReal - point(1) - point(2)- point(3), point(2), &
temp(:,16) = [point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1)] point(3), point(2), point(1), 1.0_pReal - point(1) - point(2)- point(3), &
temp(:,17) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2)] point(3), point(2), 1.0_pReal - point(1) - point(2)- point(3), point(1), &
temp(:,18) = [point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1)] point(3), 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), &
temp(:,19) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3)] point(3), 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), &
temp(:,20) = [1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2)] 1.0_pReal - point(1) - point(2)- point(3), point(1), point(2), point(3), &
temp(:,21) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3)] 1.0_pReal - point(1) - point(2)- point(3), point(1), point(3), point(2), &
temp(:,22) = [1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1)] 1.0_pReal - point(1) - point(2)- point(3), point(2), point(1), point(3), &
temp(:,23) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2)] 1.0_pReal - point(1) - point(2)- point(3), point(2), point(3), point(1), &
temp(:,24) = [1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)] 1.0_pReal - point(1) - point(2)- point(3), point(3), point(1), point(2), &
1.0_pReal - point(1) - point(2)- point(3), point(3), point(2), point(1)],[4,24])),.true.)
qPt = reshape(matmul(tetrahedron, temp),[72])
end function permutationStar1111 end function permutationStar1111
!--------------------------------------------------------------------------------------------------
!> @brief Check correctness of quadrature weights and points.
!--------------------------------------------------------------------------------------------------
subroutine selfTest
integer :: o, d, n
real(pReal), dimension(2:3), parameter :: w = [3.0_pReal,2.0_pReal]
do d = lbound(FEM_quadrature_weights,1), ubound(FEM_quadrature_weights,1)
do o = lbound(FEM_quadrature_weights(d,:),1), ubound(FEM_quadrature_weights(d,:),1)
if (dNeq(sum(FEM_quadrature_weights(d,o)%p),1.0_pReal,5e-15_pReal)) &
error stop 'quadrature weights'
enddo
enddo
do d = lbound(FEM_quadrature_points,1), ubound(FEM_quadrature_points,1)
do o = lbound(FEM_quadrature_points(d,:),1), ubound(FEM_quadrature_points(d,:),1)
n = size(FEM_quadrature_points(d,o)%p,1)/d
if (any(dNeq(sum(reshape(FEM_quadrature_points(d,o)%p,[d,n]),2),-real(n,pReal)/w(d),1.e-14_pReal))) &
error stop 'quadrature points'
enddo
enddo
end subroutine selfTest
end module FEM_quadrature end module FEM_quadrature

35
src/mesh/FEM_utilities.f90
View File

@ -23,14 +23,8 @@ module FEM_utilities
implicit none implicit none
private private
!-------------------------------------------------------------------------------------------------- logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
logical, public :: cutBack = .false. !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill real(pReal), public, protected :: wgt !< weighting factor 1/Nelems
integer, public, parameter :: maxFields = 6
integer, public :: nActiveFields = 0
!--------------------------------------------------------------------------------------------------
! grid related information information
real(pReal), public :: wgt !< weighting factor 1/Nelems
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -49,10 +43,6 @@ module FEM_utilities
COMPONENT_MECH_Z_ID COMPONENT_MECH_Z_ID
end enum end enum
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical :: &
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
@ -63,27 +53,17 @@ module FEM_utilities
end type tSolutionState end type tSolutionState
type, public :: tComponentBC type, public :: tComponentBC
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
real(pReal), allocatable, dimension(:) :: Value real(pReal), allocatable, dimension(:) :: Value
logical, allocatable, dimension(:) :: Mask logical, allocatable, dimension(:) :: Mask
end type tComponentBC end type tComponentBC
type, public :: tFieldBC type, public :: tFieldBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer :: nComponents = 0 integer :: nComponents = 0
type(tComponentBC), allocatable :: componentBC(:) type(tComponentBC), allocatable, dimension(:) :: componentBC
end type tFieldBC end type tFieldBC
type, public :: tLoadCase
real(pReal) :: time = 0.0_pReal !< length of increment
integer :: incs = 0, & !< number of increments
outputfrequency = 1, & !< frequency of result writes
logscale = 0 !< linear/logarithmic time inc flag
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: faceID
type(tFieldBC), allocatable, dimension(:) :: fieldBC
end type tLoadCase
public :: & public :: &
FEM_utilities_init, & FEM_utilities_init, &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
@ -109,8 +89,9 @@ subroutine FEM_utilities_init
integer :: structOrder !< order of displacement shape functions integer :: structOrder !< order of displacement shape functions
character(len=*), parameter :: & character(len=*), parameter :: &
PETSCDEBUG = ' -snes_view -snes_monitor ' PETSCDEBUG = ' -snes_view -snes_monitor '
PetscErrorCode :: ierr PetscErrorCode :: ierr
logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution
print'(/,a)', ' <<<+- FEM_utilities init -+>>>' print'(/,a)', ' <<<+- FEM_utilities init -+>>>'

16
src/mesh/discretization_mesh.f90
View File

@ -40,6 +40,11 @@ module discretization_mesh
mesh_maxNips !< max number of IPs in any CP element mesh_maxNips !< max number of IPs in any CP element
!!!! BEGIN DEPRECATED !!!!! !!!! BEGIN DEPRECATED !!!!!
DM, public :: geomMesh
PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
mesh_ipVolume, & !< volume associated with IP (initially!) mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!) mesh_node0 !< node x,y,z coordinates (initially!)
@ -50,11 +55,6 @@ module discretization_mesh
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
DM, public :: geomMesh
PetscInt, dimension(:), allocatable, public, protected :: &
mesh_boundaries
public :: & public :: &
discretization_mesh_init, & discretization_mesh_init, &
mesh_FEM_build_ipVolumes, & mesh_FEM_build_ipVolumes, &
@ -71,16 +71,14 @@ subroutine discretization_mesh_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
integer, allocatable, dimension(:) :: chunkPos
integer :: dimPlex, & integer :: dimPlex, &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes, & !< total number of nodes in mesh
j, l, & j, &
debug_element, debug_ip debug_element, debug_ip
PetscSF :: sf PetscSF :: sf
DM :: globalMesh DM :: globalMesh
PetscInt :: nFaceSets PetscInt :: nFaceSets
PetscInt, pointer, dimension(:) :: pFaceSets PetscInt, pointer, dimension(:) :: pFaceSets
character(len=pStringLen), dimension(:), allocatable :: fileContent
IS :: faceSetIS IS :: faceSetIS
PetscErrorCode :: ierr PetscErrorCode :: ierr
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -88,7 +86,7 @@ subroutine discretization_mesh_init(restart)
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
integer :: integrationOrder !< order of quadrature rule required integer :: integrationOrder !< order of quadrature rule required
type(tvec) :: coords_node0 type(tvec) :: coords_node0
print'(/,a)', ' <<<+- discretization_mesh init -+>>>' print'(/,a)', ' <<<+- discretization_mesh init -+>>>'

54
src/mesh/mesh_mech_FEM.f90
View File

@ -109,7 +109,7 @@ subroutine FEM_mechanical_init(fieldBC)
character(len=*), parameter :: prefix = 'mechFE_' character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: ierr PetscErrorCode :: ierr
real(pReal), dimension(3,3) :: devNull
class(tNode), pointer :: & class(tNode), pointer :: &
num_mesh num_mesh
@ -258,6 +258,7 @@ subroutine FEM_mechanical_init(fieldBC)
call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr) call DMPlexVecSetClosure(mechanical_mesh,section,solution_local,cell,px_scal,5,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)
enddo enddo
call utilities_constitutiveResponse(0.0_pReal,devNull,.true.)
end subroutine FEM_mechanical_init end subroutine FEM_mechanical_init
@ -288,8 +289,8 @@ type(tSolutionState) function FEM_mechanical_solution( &
params%timeinc = timeinc params%timeinc = timeinc
params%fieldBC = fieldBC params%fieldBC = fieldBC
call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution) call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,ierr); CHKERRQ(ierr) ! solve mechanical_snes based on solution guess (result in solution)
call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged? call SNESGetConvergedReason(mechanical_snes,reason,ierr); CHKERRQ(ierr) ! solution converged?
terminallyIll = .false. terminallyIll = .false.
if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error if (reason < 1) then ! 0: still iterating (will not occur), negative -> convergence error
@ -397,7 +398,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call Utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
ForwardData = .false. ForwardData = .false.
@ -670,7 +671,7 @@ end subroutine FEM_mechanical_converged
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate current coordinates (FEM nodal coordinates only at the moment) !> @brief Calculate current coordinates (both nodal and ip coordinates)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_updateCoords() subroutine FEM_mechanical_updateCoords()
@ -678,21 +679,35 @@ subroutine FEM_mechanical_updateCoords()
nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes) nodeCoords_linear !< nodal coordinates (dimPlex*Nnodes)
real(pReal), pointer, dimension(:,:) :: & real(pReal), pointer, dimension(:,:) :: &
nodeCoords !< nodal coordinates (3,Nnodes) nodeCoords !< nodal coordinates (3,Nnodes)
real(pReal), pointer, dimension(:,:,:) :: &
ipCoords !< ip coordinates (3,nQuadrature,mesh_NcpElems)
integer :: &
qPt, &
comp, &
qOffset, &
nOffset
DM :: dm_local DM :: dm_local
Vec :: x_local Vec :: x_local
PetscErrorCode :: ierr PetscErrorCode :: ierr
PetscInt :: dimPlex, pStart, pEnd, p, s, e PetscInt :: pStart, pEnd, p, s, e, q, &
cellStart, cellEnd, c, n
PetscSection :: section PetscSection :: section
PetscQuadrature :: mechQuad
PetscReal, dimension(:), pointer :: basisField, basisFieldDer
PetscScalar, dimension(:), pointer :: x_scal
call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr) call SNESGetDM(mechanical_snes,dm_local,ierr); CHKERRQ(ierr)
call DMGetDS(dm_local,mechQuad,ierr); CHKERRQ(ierr)
call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetLocalSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr) call DMGetDimension(dm_local,dimPlex,ierr); CHKERRQ(ierr)
! write cell vertex displacements
call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr) call DMPlexGetDepthStratum(dm_local,0,pStart,pEnd,ierr); CHKERRQ(ierr)
allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal) allocate(nodeCoords(3,pStart:pEnd-1),source=0.0_pReal)
call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) call VecGetArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
do p=pStart, pEnd-1 do p=pStart, pEnd-1
call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr) call DMPlexGetPointLocal(dm_local, p, s, e, ierr); CHKERRQ(ierr)
nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e) nodeCoords(1:dimPlex,p)=nodeCoords_linear(s+1:e)
@ -700,6 +715,31 @@ subroutine FEM_mechanical_updateCoords()
call discretization_setNodeCoords(nodeCoords) call discretization_setNodeCoords(nodeCoords)
call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr) call VecRestoreArrayF90(x_local,nodeCoords_linear,ierr); CHKERRQ(ierr)
! write ip displacements
call DMPlexGetHeightStratum(dm_local,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
call PetscDSGetTabulation(mechQuad,0,basisField,basisFieldDer,ierr); CHKERRQ(ierr)
allocate(ipCoords(3,nQuadrature,mesh_NcpElems),source=0.0_pReal)
do c=cellStart,cellEnd-1
qOffset=0
call DMPlexVecGetClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr) !< get nodal coordinates of each element
do qPt=0,nQuadrature-1
qOffset= qPt * (size(basisField)/nQuadrature)
do comp=0,dimPlex-1 !< loop over components
nOffset=0
q = comp
do n=0,nBasis-1
ipCoords(comp+1,qPt+1,c+1)=ipCoords(comp+1,qPt+1,c+1)+&
sum(basisField(qOffset+(q*dimPlex)+1:qOffset+(q*dimPlex)+dimPlex)*&
x_scal(nOffset+1:nOffset+dimPlex))
q = q+dimPlex
nOffset = nOffset+dimPlex
enddo
enddo
enddo
call DMPlexVecRestoreClosure(dm_local,section,x_local,c,x_scal,ierr); CHKERRQ(ierr)
end do
call discretization_setIPcoords(reshape(ipCoords,[3,mesh_NcpElems*nQuadrature]))
call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMRestoreLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)
end subroutine FEM_mechanical_updateCoords end subroutine FEM_mechanical_updateCoords

30
src/phase_mechanical.f90
View File

@ -229,8 +229,8 @@ module subroutine mechanical_init(materials,phases)
allocate(phase_mechanical_F0(phases%length)) allocate(phase_mechanical_F0(phases%length))
allocate(phase_mechanical_Li(phases%length)) allocate(phase_mechanical_Li(phases%length))
allocate(phase_mechanical_Li0(phases%length)) allocate(phase_mechanical_Li0(phases%length))
allocate(phase_mechanical_Lp0(phases%length))
allocate(phase_mechanical_Lp(phases%length)) allocate(phase_mechanical_Lp(phases%length))
allocate(phase_mechanical_Lp0(phases%length))
allocate(phase_mechanical_S(phases%length)) allocate(phase_mechanical_S(phases%length))
allocate(phase_mechanical_P(phases%length)) allocate(phase_mechanical_P(phases%length))
allocate(phase_mechanical_S0(phases%length)) allocate(phase_mechanical_S0(phases%length))
@ -238,20 +238,20 @@ module subroutine mechanical_init(materials,phases)
do ph = 1, phases%length do ph = 1, phases%length
Nmembers = count(material_phaseID == ph) Nmembers = count(material_phaseID == ph)
allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fe(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fi0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Fp0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Li(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers)) allocate(phase_mechanical_Li0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_Lp0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal) allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_F(ph)%data(3,3,Nmembers))
allocate(phase_mechanical_F0(ph)%data(3,3,Nmembers))
phase => phases%get(ph) phase => phases%get(ph)
mech => phase%get('mechanical') mech => phase%get('mechanical')
@ -508,7 +508,7 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
enddo LpLoop enddo LpLoop
call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & call phase_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
S, Fi_new, ph,en) S, Fi_new, ph,en)
!* update current residuum and check for convergence of loop !* update current residuum and check for convergence of loop
atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error atol_Li = max(num%rtol_crystalliteStress * max(norm2(Liguess),norm2(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error
@ -1021,7 +1021,7 @@ module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged
subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en) subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en) subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
subState0 = plasticState(ph)%State0(:,en) allocate(subState0,source=plasticState(ph)%State0(:,en))
subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en) subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en) subF0 = phase_mechanical_F0(ph)%data(1:3,1:3,en)
@ -1144,12 +1144,12 @@ module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
en = material_phaseEntry(co,ce) en = material_phaseEntry(co,ce)
call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
phase_mechanical_Fe(ph)%data(1:3,1:3,en), & phase_mechanical_Fe(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en) phase_mechanical_Fi(ph)%data(1:3,1:3,en),ph,en)
call phase_LiAndItsTangents(devNull,dLidS,dLidFi, & call phase_LiAndItsTangents(devNull,dLidS,dLidFi, &
phase_mechanical_S(ph)%data(1:3,1:3,en), & phase_mechanical_S(ph)%data(1:3,1:3,en), &
phase_mechanical_Fi(ph)%data(1:3,1:3,en), & phase_mechanical_Fi(ph)%data(1:3,1:3,en), &
ph,en) ph,en)
invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en)) invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en)) invFi = math_inv33(phase_mechanical_Fi(ph)%data(1:3,1:3,en))

77
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -18,7 +18,7 @@ submodule(phase:plastic) dislotungsten
Q_cl = 1.0_pReal !< activation energy for dislocation climb Q_cl = 1.0_pReal !< activation energy for dislocation climb
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
b_sl, & !< magnitude of Burgers vector [m] b_sl, & !< magnitude of Burgers vector [m]
D_a, & d_caron, & !< distance of spontaneous annhihilation
i_sl, & !< Adj. parameter for distance between 2 forest dislocations i_sl, & !< Adj. parameter for distance between 2 forest dislocations
f_at, & !< factor to calculate atomic volume f_at, & !< factor to calculate atomic volume
tau_Peierls, & !< Peierls stress tau_Peierls, & !< Peierls stress
@ -56,7 +56,7 @@ submodule(phase:plastic) dislotungsten
type :: tDisloTungstendependentState type :: tDisloTungstendependentState
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
Lambda_sl, & Lambda_sl, &
threshold_stress tau_pass
end type tDisloTungstendependentState end type tDisloTungstendependentState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -172,7 +172,6 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%D_0 = pl%get_asFloat('D_0') prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl') prm%Q_cl = pl%get_asFloat('Q_cl')
prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
prm%D_a = pl%get_asFloat('D_a') * prm%b_sl
prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.) prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.)
@ -191,7 +190,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
prm%B = math_expand(prm%B, N_sl) prm%B = math_expand(prm%B, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%f_at = math_expand(prm%f_at, N_sl) prm%f_at = math_expand(prm%f_at, N_sl)
prm%D_a = math_expand(prm%D_a, N_sl) prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl
! sanity checks ! 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'
@ -202,12 +201,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' 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%Q_s <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_s'
if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls' if (any(prm%tau_Peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_Peierls'
if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' D_a or b_sl' if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)'
if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl' if (any(prm%f_at <= 0.0_pReal)) extmsg = trim(extmsg)//' f_at or b_sl'
else slipActive else slipActive
rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%f_at,prm%tau_Peierls, & 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%v_0,prm%p,prm%q,prm%B,prm%h,prm%w,prm%omega, &
source = emptyRealArray) source = emptyRealArray)
allocate(prm%forestProjection(0,0)) allocate(prm%forestProjection(0,0))
@ -246,8 +246,8 @@ module function plastic_dislotungsten_init() result(myPlasticity)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,Nmembers), source=0.0_pReal) allocate(dst%tau_pass(prm%sum_N_sl,Nmembers), source=0.0_pReal)
end associate end associate
@ -316,8 +316,6 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
ph, & ph, &
en en
real(pReal) :: &
VacancyDiffusion
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_gamma_pos, dot_gamma_neg,& dot_gamma_pos, dot_gamma_neg,&
tau_pos,& tau_pos,&
@ -325,38 +323,36 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
v_cl, & v_cl, &
dot_rho_dip_formation, & dot_rho_dip_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
dip_distance d_hat
associate(prm => param(ph), stt => state(ph),& associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
dot => dotState(ph), dst => dependentState(ph))
call kinetics(Mp,T,ph,en,& call kinetics(Mp,T,ph,en,&
dot_gamma_pos,dot_gamma_neg, & dot_gamma_pos,dot_gamma_neg, &
tau_pos_out = tau_pos,tau_neg_out = tau_neg) tau_pos_out = tau_pos,tau_neg_out = tau_neg)
dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg) ! ToDo: needs to be abs dot%gamma_sl(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
where(dEq0(tau_pos)) ! ToDo: use avg of pos and neg where(dEq0(tau_pos)) ! ToDo: use avg of +/-
dot_rho_dip_formation = 0.0_pReal dot_rho_dip_formation = 0.0_pReal
dot_rho_dip_climb = 0.0_pReal dot_rho_dip_climb = 0.0_pReal
else where else where
dip_distance = math_clip(3.0_pReal*prm%mu*prm%b_sl/(16.0_pReal*PI*abs(tau_pos)), & 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_a, & ! lower limit prm%d_caron, & ! lower limit
dst%Lambda_sl(:,en)) ! upper limit dst%Lambda_sl(:,en)) ! upper limit
dot_rho_dip_formation = merge(2.0_pReal*dip_distance* stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en))/prm%b_sl, & ! ToDo: ignore region of spontaneous annihilation 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, & 0.0_pReal, &
prm%dipoleformation) prm%dipoleformation)
v_cl = (3.0_pReal*prm%mu*VacancyDiffusion*prm%f_at/(2.0_pReal*pi*kB*T)) & 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/(dip_distance+prm%D_a)) * (1.0_pReal/(d_hat+prm%d_caron))
dot_rho_dip_climb = (4.0_pReal*v_cl*stt%rho_dip(:,en))/(dip_distance-prm%D_a) ! ToDo: Discuss with Franz: Stress dependency? 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 end where
dot%rho_mob(:,en) = abs(dot%gamma_sl(:,en))/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication dot%rho_mob(:,en) = dot%gamma_sl(:,en)/(prm%b_sl*dst%Lambda_sl(:,en)) & ! multiplication
- dot_rho_dip_formation & - dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_mob(:,en)*abs(dot%gamma_sl(:,en)) ! Spontaneous annihilation of 2 single edge dislocations - (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 & dot%rho_dip(:,en) = dot_rho_dip_formation &
- (2.0_pReal*prm%D_a)/prm%b_sl*stt%rho_dip(:,en)*abs(dot%gamma_sl(:,en)) & ! Spontaneous annihilation of a single edge dislocation with a dipole constituent - (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_dip_climb
end associate end associate
@ -377,11 +373,11 @@ module subroutine dislotungsten_dependentState(ph,en)
dislocationSpacing dislocationSpacing
associate(prm => param(ph), stt => state(ph),dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
dst%threshold_stress(:,en) = prm%mu*prm%b_sl & dst%tau_pass(:,en) = prm%mu*prm%b_sl &
* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) * 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) dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
@ -416,7 +412,7 @@ module subroutine plastic_dislotungsten_results(ph,group)
if(prm%sum_N_sl>0) call results_writeDataset(dst%Lambda_sl,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(dst%Lambda_sl,group,trim(prm%output(o)), &
'mean free path for slip','m') 'mean free path for slip','m')
case('tau_pass') case('tau_pass')
if(prm%sum_N_sl>0) call results_writeDataset(dst%threshold_stress,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(dst%tau_pass,group,trim(prm%output(o)), &
'threshold stress for slip','Pa') 'threshold stress for slip','Pa')
end select end select
enddo outputsLoop enddo outputsLoop
@ -456,8 +452,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
StressRatio_p,StressRatio_pminus1, & StressRatio_p,StressRatio_pminus1, &
dvel, vel, & dvel, vel, &
tau_pos,tau_neg, & tau_pos,tau_neg, &
t_n, t_k, dtk,dtn, & t_n, t_k, dtk,dtn
needsGoodName ! ToDo: @Karo: any idea?
integer :: j integer :: j
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
@ -475,13 +470,12 @@ pure subroutine kinetics(Mp,T,ph,en, &
dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, & dot_gamma_0 => stt%rho_mob(:,en)*prm%b_sl*prm%v_0, &
effectiveLength => dst%Lambda_sl(:,en) - prm%w) effectiveLength => dst%Lambda_sl(:,en) - prm%w)
significantPositiveTau: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) significantPositiveTau: where(abs(tau_pos)-dst%tau_pass(:,en) > tol_math_check)
StressRatio = (abs(tau_pos)-dst%threshold_stress(:,en))/prm%tau_Peierls StressRatio = (abs(tau_pos)-dst%tau_pass(:,en))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) 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_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)
@ -492,7 +486,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
end where significantPositiveTau end where significantPositiveTau
if (present(ddot_gamma_dtau_pos)) then if (present(ddot_gamma_dtau_pos)) then
significantPositiveTau2: where(abs(tau_pos)-dst%threshold_stress(:,en) > tol_math_check) 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) & 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)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_pos dtk = -1.0_pReal * t_k / tau_pos
@ -505,13 +499,12 @@ pure subroutine kinetics(Mp,T,ph,en, &
end where significantPositiveTau2 end where significantPositiveTau2
endif endif
significantNegativeTau: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) significantNegativeTau: where(abs(tau_neg)-dst%tau_pass(:,en) > tol_math_check)
StressRatio = (abs(tau_neg)-dst%threshold_stress(:,en))/prm%tau_Peierls StressRatio = (abs(tau_neg)-dst%tau_pass(:,en))/prm%tau_Peierls
StressRatio_p = StressRatio** prm%p StressRatio_p = StressRatio** prm%p
StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal)
needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q)
t_n = prm%b_sl/(needsGoodName*prm%omega*effectiveLength) 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_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)
@ -522,7 +515,7 @@ pure subroutine kinetics(Mp,T,ph,en, &
end where significantNegativeTau end where significantNegativeTau
if (present(ddot_gamma_dtau_neg)) then if (present(ddot_gamma_dtau_neg)) then
significantNegativeTau2: where(abs(tau_neg)-dst%threshold_stress(:,en) > tol_math_check) 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) & 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)**(prm%p - 1.0_pReal) / prm%tau_Peierls
dtk = -1.0_pReal * t_k / tau_neg dtk = -1.0_pReal * t_k / tau_neg

286
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -16,13 +16,11 @@ submodule(phase:plastic) dislotwin
real(pReal) :: & real(pReal) :: &
mu = 1.0_pReal, & !< equivalent shear modulus mu = 1.0_pReal, & !< equivalent shear modulus
nu = 1.0_pReal, & !< equivalent shear Poisson's ratio nu = 1.0_pReal, & !< equivalent shear Poisson's ratio
D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient
Q_cl = 1.0_pReal, & !< activation energy for dislocation climb Q_cl = 1.0_pReal, & !< activation energy for dislocation climb
omega = 1.0_pReal, & !< frequency factor for dislocation climb omega = 1.0_pReal, & !< frequency factor for dislocation climb
D = 1.0_pReal, & !< grain size D = 1.0_pReal, & !< grain size
p_sb = 1.0_pReal, & !< p-exponent in shear band velocity p_sb = 1.0_pReal, & !< p-exponent in shear band velocity
q_sb = 1.0_pReal, & !< q-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity
D_a = 1.0_pReal, & !< adjustment parameter to calculate minimum dipole distance
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors L_tw = 1.0_pReal, & !< Length of twin nuclei in Burgers vectors
L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors L_tr = 1.0_pReal, & !< Length of trans nuclei in Burgers vectors
@ -42,7 +40,7 @@ submodule(phase:plastic) dislotwin
b_sl, & !< absolute length of Burgers vector [m] for each slip system b_sl, & !< absolute length of Burgers vector [m] for each slip system
b_tw, & !< absolute length of Burgers vector [m] for each twin system b_tw, & !< absolute length of Burgers vector [m] for each twin system
b_tr, & !< absolute length of Burgers vector [m] for each transformation system b_tr, & !< absolute length of Burgers vector [m] for each transformation system
Q_s,& !< activation energy for glide [J] for each slip system Q_sl,& !< activation energy for glide [J] for each slip system
v_0, & !< dislocation velocity prefactor [m/s] for each slip system v_0, & !< dislocation velocity prefactor [m/s] for each slip system
dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system
dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system
@ -55,7 +53,8 @@ submodule(phase:plastic) dislotwin
s, & !< s-exponent in trans nucleation rate s, & !< s-exponent in trans nucleation rate
tau_0, & !< strength due to elements in solid solution tau_0, & !< strength due to elements in solid solution
gamma_char, & !< characteristic shear for twins gamma_char, & !< characteristic shear for twins
B !< drag coefficient B, & !< drag coefficient
d_caron !< distance of spontaneous annhihilation
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
h_sl_sl, & !< components of slip-slip interaction matrix h_sl_sl, & !< components of slip-slip interaction matrix
h_sl_tw, & !< components of slip-twin interaction matrix h_sl_tw, & !< components of slip-twin interaction matrix
@ -206,7 +205,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
rho_dip_0 = pl%get_as1dFloat('rho_dip_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%v_0 = pl%get_as1dFloat('v_0', requiredSize=size(N_sl))
prm%b_sl = pl%get_as1dFloat('b_sl', 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)) prm%Q_sl = pl%get_as1dFloat('Q_sl', requiredSize=size(N_sl))
prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl)) prm%i_sl = pl%get_as1dFloat('i_sl', requiredSize=size(N_sl))
prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl)) prm%p = pl%get_as1dFloat('p_sl', requiredSize=size(N_sl))
prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl)) prm%q = pl%get_as1dFloat('q_sl', requiredSize=size(N_sl))
@ -214,9 +213,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), & prm%B = pl%get_as1dFloat('B', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal, i=1,size(N_sl))]) defaultVal=[(0.0_pReal, i=1,size(N_sl))])
prm%D_a = pl%get_asFloat('D_a') prm%Q_cl = pl%get_asFloat('Q_cl')
prm%D_0 = pl%get_asFloat('D_0')
prm%Q_cl = pl%get_asFloat('Q_cl')
prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.)
prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.)
@ -231,28 +228,29 @@ module function plastic_dislotwin_init() result(myPlasticity)
rho_dip_0 = math_expand(rho_dip_0, N_sl) rho_dip_0 = math_expand(rho_dip_0, N_sl)
prm%v_0 = math_expand(prm%v_0, N_sl) prm%v_0 = math_expand(prm%v_0, N_sl)
prm%b_sl = math_expand(prm%b_sl, N_sl) prm%b_sl = math_expand(prm%b_sl, N_sl)
prm%Q_s = math_expand(prm%Q_s, N_sl) prm%Q_sl = math_expand(prm%Q_sl, N_sl)
prm%i_sl = math_expand(prm%i_sl, N_sl) prm%i_sl = math_expand(prm%i_sl, N_sl)
prm%p = math_expand(prm%p, N_sl) prm%p = math_expand(prm%p, N_sl)
prm%q = math_expand(prm%q, N_sl) prm%q = math_expand(prm%q, N_sl)
prm%tau_0 = math_expand(prm%tau_0, N_sl) prm%tau_0 = math_expand(prm%tau_0, N_sl)
prm%B = math_expand(prm%B, N_sl) prm%B = math_expand(prm%B, N_sl)
prm%d_caron = pl%get_asFloat('D_a') * prm%b_sl
! sanity checks ! sanity checks
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 ( 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_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(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%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%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%Q_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' Q_sl'
if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl' if (any(prm%i_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' i_sl'
if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B'
if (any(prm%d_caron < 0.0_pReal)) extmsg = trim(extmsg)//' d_caron(D_a,b_sl)'
if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p_sl'
if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q_sl'
else slipActive else slipActive
rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray rho_mob_0 = emptyRealArray; rho_dip_0 = emptyRealArray
allocate(prm%b_sl,prm%Q_s,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray) allocate(prm%b_sl,prm%Q_sl,prm%v_0,prm%i_sl,prm%p,prm%q,prm%B,source=emptyRealArray)
allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0)) allocate(prm%forestProjection(0,0),prm%h_sl_sl(0,0))
endif slipActive endif slipActive
@ -516,7 +514,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
integer :: i,k,l,m,n integer :: i,k,l,m,n
real(pReal) :: & real(pReal) :: &
f_unrotated,StressRatio_p,& f_unrotated,StressRatio_p,&
BoltzmannRatio, & E_kB_T, &
ddot_gamma_dtau, & ddot_gamma_dtau, &
tau tau
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
@ -586,7 +584,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
shearBandingContribution: if(dNeq0(prm%v_sb)) then shearBandingContribution: if(dNeq0(prm%v_sb)) then
BoltzmannRatio = prm%E_sb/(kB*T) E_kB_T = prm%E_sb/(kB*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
do i = 1,6 do i = 1,6
@ -596,8 +594,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
significantShearBandStress: if (abs(tau) > tol_math_check) then significantShearBandStress: if (abs(tau) > tol_math_check) then
StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) dot_gamma_sb = sign(prm%v_sb*exp(-E_kB_T*(1-StressRatio_p)**prm%q_sb), tau)
ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb & ddot_gamma_dtau = abs(dot_gamma_sb)*E_kB_T*prm%p_sb*prm%q_sb/prm%xi_sb &
* (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) & * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
* (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
@ -631,7 +629,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
integer :: i integer :: i
real(pReal) :: & real(pReal) :: &
f_unrotated, & f_unrotated, &
rho_dip_distance, & d_hat, &
v_cl, & !< climb velocity v_cl, & !< climb velocity
tau, & tau, &
sigma_cl, & !< climb stress sigma_cl, & !< climb stress
@ -639,70 +637,67 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
dot_rho_dip_formation, & dot_rho_dip_formation, &
dot_rho_dip_climb, & dot_rho_dip_climb, &
rho_dip_distance_min, &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
dot_gamma_tw dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
associate(prm => param(ph), stt => state(ph), &
dot => dotState(ph), dst => dependentState(ph))
f_unrotated = 1.0_pReal & associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
- sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en))
call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) f_unrotated = 1.0_pReal &
dot%gamma_sl(:,en) = abs(dot_gamma_sl) - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en))
rho_dip_distance_min = prm%D_a*prm%b_sl call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
dot%gamma_sl(:,en) = abs(dot_gamma_sl)
slipState: do i = 1, prm%sum_N_sl slipState: do i = 1, prm%sum_N_sl
tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then significantSlipStress: if (dEq0(tau) .or. prm%omitDipoles) then
dot_rho_dip_formation(i) = 0.0_pReal dot_rho_dip_formation(i) = 0.0_pReal
dot_rho_dip_climb(i) = 0.0_pReal
else significantSlipStress
rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,en))
rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i))
dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) &
* stt%rho_mob(i,en)*abs(dot_gamma_sl(i))
if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else else significantSlipStress
! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 d_hat = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau))
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) d_hat = math_clip(d_hat, right = dst%Lambda_sl(i,en))
b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) & d_hat = math_clip(d_hat, left = prm%d_caron(i))
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
1.0_pReal, &
prm%ExtendedDislocations)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) & dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) &
/ (rho_dip_distance-rho_dip_distance_min(i)) * stt%rho_mob(i,en)*abs(dot_gamma_sl(i))
endif
endif significantSlipStress
enddo slipState
dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & if (dEq(d_hat,prm%d_caron(i))) then
- dot_rho_dip_formation & dot_rho_dip_climb(i) = 0.0_pReal
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) else
! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
b_d = merge(24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu) &
* (prm%Gamma_sf(1) + prm%Gamma_sf(2) * T) / (prm%mu*prm%b_sl(i)), &
1.0_pReal, &
prm%ExtendedDislocations)
v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Q_cl/(kB*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal)
dot%rho_dip(:,en) = dot_rho_dip_formation & dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,en) &
- 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & / (d_hat-prm%d_caron(i))
- dot_rho_dip_climb endif
endif significantSlipStress
enddo slipState
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) &
dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char - dot_rho_dip_formation &
- 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl)
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) dot%rho_dip(:,en) = dot_rho_dip_formation &
dot%f_tr(:,en) = f_unrotated*dot_gamma_tr - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) &
- dot_rho_dip_climb
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
end associate end associate
@ -763,19 +758,17 @@ module subroutine dislotwin_dependentState(T,ph,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%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
!* threshold stress for growing twin/martensite !* threshold stress for growing twin/martensite
if(prm%sum_N_tw == prm%sum_N_sl) & dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) & + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_tw)
+ 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct? dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
if(prm%sum_N_tr == prm%sum_N_sl) & + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_tr) &
dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) & + prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
+ 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+ prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0) dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the Burgers vector for slip
@ -867,7 +860,7 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
tau, & tau, &
stressRatio, & stressRatio, &
StressRatio_p, & StressRatio_p, &
BoltzmannRatio, & Q_kB_T, &
v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned)
v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned)
dV_wait_inverse_dTau, & dV_wait_inverse_dTau, &
@ -876,33 +869,34 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
tau_eff !< effective resolved stress tau_eff !< effective resolved stress
integer :: i integer :: i
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)] tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)]
tau_eff = abs(tau)-dst%tau_pass(:,en) tau_eff = abs(tau)-dst%tau_pass(:,en)
significantStress: where(tau_eff > tol_math_check) significantStress: where(tau_eff > tol_math_check)
stressRatio = tau_eff/prm%tau_0 stressRatio = tau_eff/prm%tau_0
StressRatio_p = stressRatio** prm%p StressRatio_p = stressRatio** prm%p
BoltzmannRatio = prm%Q_s/(kB*T) Q_kB_T = prm%Q_sl/(kB*T)
v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_wait_inverse = prm%v_0**(-1.0_pReal) * exp(Q_kB_T*(1.0_pReal-StressRatio_p)** prm%q)
v_run_inverse = prm%B/(tau_eff*prm%b_sl) 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) dot_gamma_sl = sign(stt%rho_mob(:,en)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau)
dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * Q_kB_T &
* (stressRatio**(prm%p-1.0_pReal)) & * (stressRatio**(prm%p-1.0_pReal)) &
* (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) &
/ prm%tau_0 / prm%tau_0
dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff
dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) &
/ (v_wait_inverse+v_run_inverse)**2.0_pReal / (v_wait_inverse+v_run_inverse)**2.0_pReal
ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,en)*prm%b_sl
else where significantStress else where significantStress
dot_gamma_sl = 0.0_pReal dot_gamma_sl = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate
@ -945,34 +939,35 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
integer :: i,s1,s2 integer :: i,s1,s2
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_tw do i = 1, prm%sum_N_tw
tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) tau(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
isFCC: if (prm%fccTwinTransNucleation) then isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i) s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i) s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct? if (tau(i) < dst%tau_r_tw(i,en)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
(prm%L_tw*prm%b_sl(i))*& (prm%L_tw*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i)))) ! P_ncs (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,en)-tau(i))))
else else
Ndot0=0.0_pReal Ndot0=0.0_pReal
end if end if
else isFCC else isFCC
Ndot0=prm%dot_N_0_tw(i) Ndot0=prm%dot_N_0_tw(i)
endif isFCC endif isFCC
enddo enddo
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r StressRatio_r = (dst%tau_hat_tw(:,en)/tau)**prm%r
dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r) dot_gamma_tw = prm%gamma_char * dst%V_tw(:,en) * Ndot0*exp(-StressRatio_r)
ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r ddot_gamma_dtau = (dot_gamma_tw*prm%r/tau)*StressRatio_r
else where significantStress else where significantStress
dot_gamma_tw = 0.0_pReal dot_gamma_tw = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate
@ -1011,36 +1006,37 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
Ndot0, & Ndot0, &
stressRatio_s, & stressRatio_s, &
ddot_gamma_dtau ddot_gamma_dtau
integer :: i,s1,s2 integer :: i,s1,s2
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
do i = 1, prm%sum_N_tr do i = 1, prm%sum_N_tr
tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) tau(i) = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
isFCC: if (prm%fccTwinTransNucleation) then isFCC: if (prm%fccTwinTransNucleation) then
s1=prm%fcc_twinNucleationSlipPair(1,i) s1=prm%fcc_twinNucleationSlipPair(1,i)
s2=prm%fcc_twinNucleationSlipPair(2,i) s2=prm%fcc_twinNucleationSlipPair(2,i)
if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct? if (tau(i) < dst%tau_r_tr(i,en)) then ! ToDo: correct?
Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+& Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,en)+stt%rho_dip(s2,en))+&
abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/& ! ToDo: MD: it would be more consistent to use shearrates from state abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,en)+stt%rho_dip(s1,en)))/&
(prm%L_tr*prm%b_sl(i))*& (prm%L_tr*prm%b_sl(i))*&
(1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i)))) ! P_ncs (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,en)-tau(i))))
else else
Ndot0=0.0_pReal Ndot0=0.0_pReal
end if end if
else isFCC else isFCC
Ndot0=prm%dot_N_0_tr(i) Ndot0=prm%dot_N_0_tr(i)
endif isFCC endif isFCC
enddo enddo
significantStress: where(tau > tol_math_check) significantStress: where(tau > tol_math_check)
StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s StressRatio_s = (dst%tau_hat_tr(:,en)/tau)**prm%s
dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s) dot_gamma_tr = dst%V_tr(:,en) * Ndot0*exp(-StressRatio_s)
ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s
else where significantStress else where significantStress
dot_gamma_tr = 0.0_pReal dot_gamma_tr = 0.0_pReal
ddot_gamma_dtau = 0.0_pReal ddot_gamma_dtau = 0.0_pReal
end where significantStress end where significantStress
end associate end associate

4
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -359,11 +359,11 @@ module subroutine plastic_kinehardening_results(ph,group)
case ('xi') case ('xi')
if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), &
'resistance against plastic slip','Pa') 'resistance against plastic slip','Pa')
case ('tau_b') !ToDo: chi case ('chi')
if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), &
'back stress','Pa') 'back stress','Pa')
case ('sgn(gamma)') case ('sgn(gamma)')
if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int if(prm%sum_N_sl>0) call results_writeDataset(int(stt%sgn_gamma),group,trim(prm%output(o)), &
'sense of shear','1') 'sense of shear','1')
case ('chi_0') case ('chi_0')
if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), & if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), &

5
src/prec.f90
View File

@ -27,11 +27,6 @@ module prec
real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) real(pReal), parameter :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
type :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p
end type group_float
type :: tState type :: tState
integer :: & integer :: &
sizeState = 0, & !< size of state sizeState = 0, & !< size of state

5
src/quit.f90
View File

@ -6,7 +6,10 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine quit(stop_id) subroutine quit(stop_id)
#include <petsc/finclude/petscsys.h> #include <petsc/finclude/petscsys.h>
use PetscSys use PETScSys
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>14) && !defined(PETSC_HAVE_MPI_F90MODULE_VISIBILITY)
use MPI_f08
#endif
use HDF5 use HDF5
implicit none implicit none