From 482ba98b18fc92e8bec7a63b55ab1cc1027e3821 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 2 Aug 2023 13:27:59 +0200 Subject: [PATCH] mesh numerics changes --- PRIVATE | 2 +- examples/config/numerics.yaml | 21 ++++++++++----------- src/mesh/DAMASK_mesh.f90 | 16 +++++++++------- src/mesh/FEM_utilities.f90 | 20 +++++++++++--------- src/mesh/discretization_mesh.f90 | 4 +++- src/mesh/mesh_mech_FEM.f90 | 20 +++++++++----------- 6 files changed, 43 insertions(+), 40 deletions(-) diff --git a/PRIVATE b/PRIVATE index 8c05965ef..9a067cb27 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8c05965ef4437598898f467a213ffa88938e860a +Subproject commit 9a067cb2757df91b9242414779a8b52ac9e5ee02 diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index 487b6cf88..247ae21c5 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -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: diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 4bfd791f5..133a98997 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -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 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index d30b223f0..fb54581f1 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -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) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 5cd12549e..93476b0fb 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -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) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 43a1b0579..27ff89f9c 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -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