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

added phi_min back to example; more consistency in code syntax

This commit is contained in:
Philip Eisenlohr 2023-09-06 15:11:48 -04:00
parent 925da8aad4
commit af955df891
5 changed files with 15 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
examples/config/numerics.yaml
View File

@ -10,6 +10,7 @@ solver:
N_iter_max: 100 # maximum iteration number
eps_abs_phi: 1.0e-2 # absolute tolerance for damage evolution
eps_rel_phi: 1.0e-6 # relative tolerance for damage evolution
phi_min: 1.0e-6 # residual integrity
thermal:
N_iter_max: 100 # maximum iteration number
@ -42,7 +43,7 @@ solver:
N_staggered_iter_max: 10 # max number of field level staggered iterations
p_s: 2 # order of displacement shape functions
p_i: 2 # order of quadrature rule required
bbarstabilisation: false
bbarstabilization: false
mechanical:
N_iter_max: 250 # Maximum iteration number

4
src/Marc/discretization_Marc.f90
View File

@ -76,9 +76,9 @@ subroutine discretization_Marc_init
print'(/,a)', ' <<<+- discretization_Marc init -+>>>'; flush(6)
num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
num_commercialFEM => num_solver%get_dict('Marc', defaultVal = emptyDict)
num_commercialFEM => num_solver%get_dict('Marc',defaultVal=emptyDict)
mesh_unitlength = num_commercialFEM%get_asReal('unit_length',defaultVal=1.0_pREAL) ! set physical extent of a length unit in mesh
if (mesh_unitlength <= 0.0_pREAL) call IO_error(301,'unitlength')
if (mesh_unitlength <= 0.0_pREAL) call IO_error(301,'unit_length')
call inputRead(elem,node0_elem,connectivity_elem,materialAt)
nElems = size(connectivity_elem,2)

4
src/mesh/DAMASK_mesh.f90
View File

@ -91,8 +91,8 @@ program DAMASK_mesh
!---------------------------------------------------------------------
! reading field information from numerics file and do sanity checks
num_solver => config_numerics%get_dict('solver', defaultVal=emptyDict)
num_mesh => num_solver%get_dict('mesh', defaultVal=emptyDict)
num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
num_mesh => num_solver%get_dict('mesh',defaultVal=emptyDict)
stagItMax = num_mesh%get_asInt('N_staggered_iter_max',defaultVal=10)
maxCutBack = num_mesh%get_asInt('N_cutback_max',defaultVal=3)

2
src/mesh/discretization_mesh.f90
View File

@ -102,7 +102,7 @@ subroutine discretization_mesh_init(restart)
! read numerics parameter
num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
num_mesh => num_solver%get_dict('mesh',defaultVal=emptyDict)
p_i = num_mesh%get_asInt('p_i',defaultVal = 2)
p_i = num_mesh%get_asInt('p_i',defaultVal=2)
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16)
call DMPlexCreateFromFile(PETSC_COMM_WORLD,CLI_geomFile,'n/a',PETSC_TRUE,globalMesh,err_PETSc)

16
src/mesh/mesh_mech_FEM.f90
View File

@ -47,7 +47,7 @@ module mesh_mechanical_FEM
p_i, & !< integration order (quadrature rule)
itmax
logical :: &
BBarStabilisation
BBarStabilization
real(pREAL) :: &
eps_struct_atol, & !< absolute tolerance for mechanical equilibrium
eps_struct_rtol !< relative tolerance for mechanical equilibrium
@ -135,12 +135,12 @@ subroutine FEM_mechanical_init(fieldBC,num_mesh)
! read numerical parametes and do sanity checks
num_mech => num_mesh%get_dict('mechanical', defaultVal=emptyDict)
num%p_i = int(num_mesh%get_asInt('p_i',defaultVal = 2),pPETSCINT)
num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
num%p_i = int(num_mesh%get_asInt('p_i',defaultVal=2),pPETSCINT)
num%BBarStabilization = num_mesh%get_asBool('bbarstabilization',defaultVal=.false.)
num%itmax = int(num_mech%get_asInt('N_iter_max',defaultVal=250),pPETSCINT)
num%eps_struct_atol = num_mech%get_asReal('eps_abs_div(P)', defaultVal = 1.0e-10_pREAL)
num%eps_struct_rtol = num_mech%get_asReal('eps_rel_div(P)', defaultVal = 1.0e-4_pREAL)
num%eps_struct_atol = num_mech%get_asReal('eps_abs_div(P)', defaultVal=1.0e-10_pREAL)
num%eps_struct_rtol = num_mech%get_asReal('eps_rel_div(P)', defaultVal=1.0e-4_pREAL)
if (num%itmax <= 1) call IO_error(301,ext_msg='N_iter_max')
if (num%eps_struct_rtol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_rel_div(P)')
@ -439,7 +439,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
end do
homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
end do
if (num%BBarStabilisation) then
if (num%BBarStabilization) then
detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature,pREAL))
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
@ -590,7 +590,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,m), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1_pPETSCINT)
if (num%BBarStabilisation) then
if (num%BBarStabilization) then
F(1:dimPlex,1:dimPlex) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex])
FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0_pREAL/real(dimPlex,pREAL))
@ -606,7 +606,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
K_eA = K_eA + matmul(transpose(BMat),MatA)
end if
end do
if (num%BBarStabilisation) then
if (num%BBarStabilization) then
FInv = math_inv33(FAvg)
K_e = K_eA*math_det33(FAvg/real(nQuadrature,pREAL))**(1.0_pREAL/real(dimPlex,pREAL)) + &
(matmul(matmul(transpose(BMatAvg), &