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

mesh numerics changes

This commit is contained in:
Sharan Roongta 2023-08-02 13:27:59 +02:00
parent 485839b651
commit 482ba98b18
6 changed files with 43 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 8c05965ef4437598898f467a213ffa88938e860a
Subproject commit 9a067cb2757df91b9242414779a8b52ac9e5ee02

21
examples/config/numerics.yaml
View File

@ -40,7 +40,7 @@ solver:
eps_rel_P: 1.0e-3 # relative tolerance for fulfillment of stress BC
N_iter_min: 1 # minimum iteration number
N_iter_max: 100 # maximum iteration number
update_gamma: false # Update Gamma-operator with current dPdF (not possible if memory_efficient=1)
update_gamma: false # Update Gamma-operator with current dPdF (not possible if memory_efficient=1)
FFT:
memory_efficient: True # Precalculate Gamma-operator (81 double per point)
@ -54,16 +54,15 @@ solver:
eps_abs_curl(F): 1.0e-10 # absolute tolerance for fulfillment of strain compatibility
eps_rel_curl(F): 5.0e-4 # relative tolerance for fulfillment of strain compatibility
mesh:
maxCutBack: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
maxStaggeredIter: 10 # max number of field level staggered iterations
structorder: 2 # order of displacement shape functions (when mesh is defined)
bbarstabilisation: false
integrationorder: 2 # order of quadrature rule required (when mesh is defined)
itmax: 250 # Maximum iteration number
itmin: 2 # Minimum iteration number
eps_struct_atol: 1.0e-10 # absolute tolerance for mechanical equilibrium
eps_struct_rtol: 1.0e-4 # relative tolerance for mechanical equilibrium
mesh:
N_cutback_max: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc)
N_staggered_iter_max: 10 # max number of field level staggered iterations
p_s: 2 # order of displacement shape functions (when mesh is defined)
bbarstabilisation: false
p_i: 2 # order of quadrature rule required (when mesh is defined)
N_iter_max: 250 # Maximum iteration number
eps_abs_div(P): 1.0e-10 # absolute tolerance for mechanical equilibrium
eps_rel_div(P): 1.0e-4 # relative tolerance for mechanical equilibrium
phase:
mechanical:

16
src/mesh/DAMASK_mesh.f90
View File

@ -66,6 +66,7 @@ program DAMASK_mesh
stagIter, &
component
type(tDict), pointer :: &
num_solver, &
num_mesh
character(len=pSTRLEN), dimension(:), allocatable :: fileContent
character(len=pSTRLEN) :: &
@ -90,12 +91,13 @@ program DAMASK_mesh
!---------------------------------------------------------------------
! reading field information from numerics file and do sanity checks
num_mesh => config_numerics%get_dict('mesh', defaultVal=emptyDict)
stagItMax = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10)
maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3)
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)
if (stagItMax < 0) call IO_error(301,ext_msg='maxStaggeredIter')
if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack')
if (stagItMax < 0) call IO_error(301,ext_msg='N_staggered_iter_max')
if (maxCutBack < 0) call IO_error(301,ext_msg='N_cutback_max')
! reading basic information from load case file and allocate data structure containing load cases
call DMGetDimension(geomMesh,dimPlex,err_PETSc) !< dimension of mesh (2D or 3D)
@ -229,8 +231,8 @@ program DAMASK_mesh
!--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers
call FEM_Utilities_init()
call FEM_mechanical_init(loadCases(1)%fieldBC(1))
call FEM_Utilities_init(num_mesh)
call FEM_mechanical_init(loadCases(1)%fieldBC(1),num_mesh)
call config_numerics_deallocate()
if (worldrank == 0) then

20
src/mesh/FEM_utilities.f90
View File

@ -90,11 +90,14 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief Allocate all neccessary fields.
!--------------------------------------------------------------------------------------------------
subroutine FEM_utilities_init
subroutine FEM_utilities_init(num_mesh)
character(len=pSTRLEN) :: petsc_optionsOrder
type(tDict), pointer :: &
type(tDict), pointer, intent(in) :: &
num_mesh
character(len=pSTRLEN) :: petsc_optionsOrder
character(len=:), allocatable :: &
extmsg, &
petsc_options
integer :: &
p_s, & !< order of shape functions
p_i !< integration order (quadrature rule)
@ -103,8 +106,6 @@ subroutine FEM_utilities_init
print'(/,1x,a)', '<<<+- FEM_utilities init -+>>>'
num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict)
p_s = num_mesh%get_asInt('p_s',defaultVal = 2)
p_i = num_mesh%get_asInt('p_i',defaultVal = p_s)
@ -117,12 +118,13 @@ subroutine FEM_utilities_init
call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc)
CHKERRQ(err_PETSc)
CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls &
petsc_options = misc_prefixOptions('-mechanical_snes_type newtonls &
&-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew &
&-mechanical_snes_ksp_ew_rtol0 0.01 -mechanical_snes_ksp_ew_rtolmax 0.01 &
&-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25', err_PETSc)
CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asStr('PETSc_options',defaultVal=''),err_PETSc)
&-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25' // &
num_mesh%get_asStr('PETSc_options',defaultVal=''), 'mechanical_')
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc)
CHKERRQ(err_PETSc)
write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),err_PETSc)

4
src/mesh/discretization_mesh.f90
View File

@ -89,6 +89,7 @@ subroutine discretization_mesh_init(restart)
PetscInt, dimension(:), allocatable :: &
materialAt
type(tDict), pointer :: &
num_solver, &
num_mesh
integer :: p_i, dim !< integration order (quadrature rule)
type(tvec) :: coords_node0
@ -99,7 +100,8 @@ subroutine discretization_mesh_init(restart)
!--------------------------------------------------------------------------------
! read numerics parameter
num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict)
num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict)
p_i = num_mesh%get_asInt('p_i',defaultVal = 2)
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16)

20
src/mesh/mesh_mech_FEM.f90
View File

@ -94,9 +94,10 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data
!--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_init(fieldBC)
subroutine FEM_mechanical_init(fieldBC,num_mesh)
type(tFieldBC), intent(in) :: fieldBC
type(tFieldBC), intent(in) :: fieldBC
type(tDict), pointer, intent(in) :: num_mesh
DM :: mechanical_mesh
PetscFE :: mechFE
@ -126,23 +127,20 @@ subroutine FEM_mechanical_init(fieldBC)
character(len=*), parameter :: prefix = 'mechFE_'
PetscErrorCode :: err_PETSc
real(pREAL), dimension(3,3) :: devNull
type(tDict), pointer :: &
num_mesh
print'(/,1x,a)', '<<<+- FEM_mech init -+>>>'; flush(IO_STDOUT)
!-----------------------------------------------------------------------------
! read numerical parametes and do sanity checks
num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict)
num%p_i = int(num_mesh%get_asInt('p_i',defaultVal = 2),pPETSCINT)
num%itmax = int(num_mesh%get_asInt('itmax',defaultVal=250),pPETSCINT)
num%itmax = int(num_mesh%get_asInt('N_iter_max',defaultVal=250),pPETSCINT)
num%BBarStabilisation = num_mesh%get_asBool('bbarstabilisation',defaultVal = .false.)
num%eps_struct_atol = num_mesh%get_asReal('eps_struct_atol', defaultVal = 1.0e-10_pREAL)
num%eps_struct_rtol = num_mesh%get_asReal('eps_struct_rtol', defaultVal = 1.0e-4_pREAL)
num%eps_struct_atol = num_mesh%get_asReal('eps_abs_div(P)', defaultVal = 1.0e-10_pREAL)
num%eps_struct_rtol = num_mesh%get_asReal('eps_rel_div(P)', defaultVal = 1.0e-4_pREAL)
if (num%itmax <= 1) call IO_error(301,ext_msg='itmax')
if (num%eps_struct_rtol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_struct_rtol')
if (num%eps_struct_atol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_struct_atol')
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)')
if (num%eps_struct_atol <= 0.0_pREAL) call IO_error(301,ext_msg='eps_abs_div(P)')
!--------------------------------------------------------------------------------------------------
! Setup FEM mech mesh