From 482ba98b18fc92e8bec7a63b55ab1cc1027e3821 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 2 Aug 2023 13:27:59 +0200 Subject: [PATCH 01/10] 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 From ac63ee36008ecbe6ace8a2509649a42a271c5273 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 2 Aug 2023 14:08:25 +0200 Subject: [PATCH 02/10] marc and other minor changes --- examples/config/numerics.yaml | 105 +++++++++++++++--------------- src/Marc/discretization_Marc.f90 | 6 +- src/grid/grid_damage_spectral.f90 | 4 +- src/mesh/FEM_utilities.f90 | 1 - 4 files changed, 59 insertions(+), 57 deletions(-) diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index 247ae21c5..8e6a91389 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -1,5 +1,6 @@ # Available numerical parameters # Case sensitive keys +--- homogenization: mechanical: @@ -21,74 +22,74 @@ homogenization: solver: grid: - N_staggered_iter_max: 10 # max number of field level staggered iterations - 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 + N_cutback_max: 3 # maximum cut back level (0: 1, 1: 0.5, 2: 0.25, etc) damage: - 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 + 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 # non-zero residual damage thermal: - N_iter_max: 100 # maximum iteration number - eps_abs_T: 1.0e-2 # absolute tolerance for thermal equilibrium - eps_rel_T: 1.0e-6 # relative tolerance for thermal equilibrium + N_iter_max: 100 # maximum iteration number + eps_abs_T: 1.0e-2 # absolute tolerance for thermal equilibrium + eps_rel_T: 1.0e-6 # relative tolerance for thermal equilibrium mechanical: - eps_abs_div(P): 1.0e-4 # absolute tolerance for fulfillment of stress equilibrium - eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium - eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC - 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) + eps_abs_div(P): 1.0e-4 # absolute tolerance for fulfillment of stress equilibrium + eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium + eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC + 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) FFT: - memory_efficient: True # Precalculate Gamma-operator (81 double per point) - divergence_correction: size+grid # Use size-independent divergence criterion - derivative: continuous # Approximation used for derivatives in Fourier space - FFTW_plan_mode: FFTW_MEASURE # planing-rigor flag, see manual on www.fftw.org - FFTW_timelimit: -1.0 # timelimit of plan creation for FFTW, see manual on www.fftw.org. -1.0: disable timelimit + memory_efficient: True # Precalculate Gamma-operator (81 double per point) + divergence_correction: size+grid # Use size-independent divergence criterion + derivative: continuous # Approximation used for derivatives in Fourier space + FFTW_plan_mode: FFTW_MEASURE # planing-rigor flag, see manual on www.fftw.org + FFTW_timelimit: -1.0 # timelimit of plan creation for FFTW, see manual on www.fftw.org. -1.0: disable timelimit + alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme + beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme + 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 PETSc_options: -snes_type ngmres -snes_ngmres_anderson # PETSc solver options - alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme - beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme - 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: - 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) + 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 + 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 + + Marc: + unit_length: 1.0 # physical length of one computational length unit phase: mechanical: - r_cutback_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation - r_cutback: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) - r_increase: 1.5 # factor to increase size of next step when previous step converged in phase state calculation - eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) - N_iter_state_max: 10 # state loop limit + r_cutback_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation + r_cutback: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) + r_increase: 1.5 # factor to increase size of next step when previous step converged in phase state calculation + eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) + N_iter_state_max: 10 # state loop limit plastic: - r_linesearch_Lp: 0.5 # factor to decrease the step due to non-convergence in Lp calculation - eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum - eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum - N_iter_Lp_max: 40 # stress loop limit for Lp - f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp - integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) + r_linesearch_Lp: 0.5 # factor to decrease the step due to non-convergence in Lp calculation + eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum + eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum + N_iter_Lp_max: 40 # stress loop limit for Lp + f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp + integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) eigen: - r_linesearch_Li: 0.5 # factor to decrease the step due to non-convergence in Li calculation - eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum - eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum - N_iter_Li_max: 40 # stress loop limit for Li - f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li - -commercialFEM: - unitlength: 1 # physical length of one computational length unit + r_linesearch_Li: 0.5 # factor to decrease the step due to non-convergence in Li calculation + eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum + eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum + N_iter_Li_max: 40 # stress loop limit for Li + f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li generic: - random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed. - phi_min: 1.0e-6 # non-zero residual damage. + random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed. diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90 index 3d9778e0e..b38e3b2ac 100644 --- a/src/Marc/discretization_Marc.f90 +++ b/src/Marc/discretization_Marc.f90 @@ -69,13 +69,15 @@ subroutine discretization_Marc_init unscaledNormals type(tDict), pointer :: & + num_solver, & num_commercialFEM print'(/,a)', ' <<<+- discretization_Marc init -+>>>'; flush(6) - num_commercialFEM => config_numerics%get_dict('commercialFEM',defaultVal = emptyDict) - mesh_unitlength = num_commercialFEM%get_asReal('unitlength',defaultVal=1.0_pREAL) ! set physical extent of a length unit in mesh + num_solver => config_numerics%get_asDict('solver',defaultVal=emptyDict) + num_commercialFEM => num_solver%get_asDict('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') call inputRead(elem,node0_elem,connectivity_elem,materialAt) diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 48023cc54..648f9ebbb 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -102,7 +102,7 @@ subroutine grid_damage_spectral_init(num_grid) ! read numerical parameters and do sanity checks num_grid_damage => num_grid%get_dict('damage',defaultVal=emptyDict) - num%itmax = num_grid_damage%get_asInt('N_iter_max', defaultVal=100) + num%itmax = num_grid_damage%get_asInt ('N_iter_max', defaultVal=100) num%eps_damage_atol = num_grid_damage%get_asReal('eps_abs_phi',defaultVal=1.0e-2_pREAL) num%eps_damage_rtol = num_grid_damage%get_asReal('eps_rel_phi',defaultVal=1.0e-6_pREAL) num%phi_min = num_grid_damage%get_asReal('phi_min', defaultVal=1.0e-6_pREAL) @@ -110,7 +110,7 @@ subroutine grid_damage_spectral_init(num_grid) extmsg = '' if (num%eps_damage_atol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_abs_phi' if (num%eps_damage_rtol <= 0.0_pREAL) extmsg = trim(extmsg)//' eps_rel_phi' - if (num%phi_min < 0.0_pREAL) extmsg = trim(extmsg)//' phi_min' + if (num%phi_min <= 0.0_pREAL) extmsg = trim(extmsg)//' phi_min' if (num%itmax < 1) extmsg = trim(extmsg)//' N_iter_max' if (extmsg /= '') call IO_error(301,ext_msg=trim(extmsg)) diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index fb54581f1..d6b0c7b77 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -96,7 +96,6 @@ subroutine FEM_utilities_init(num_mesh) num_mesh character(len=pSTRLEN) :: petsc_optionsOrder character(len=:), allocatable :: & - extmsg, & petsc_options integer :: & p_s, & !< order of shape functions From c960f8ed582334f48bd75e39161adead13028e0b Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 2 Aug 2023 14:20:17 +0200 Subject: [PATCH 03/10] space needed to add new option + marc changes --- src/Marc/discretization_Marc.f90 | 4 ++-- src/mesh/FEM_utilities.f90 | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90 index b38e3b2ac..a5094b7fa 100644 --- a/src/Marc/discretization_Marc.f90 +++ b/src/Marc/discretization_Marc.f90 @@ -75,8 +75,8 @@ subroutine discretization_Marc_init print'(/,a)', ' <<<+- discretization_Marc init -+>>>'; flush(6) - num_solver => config_numerics%get_asDict('solver',defaultVal=emptyDict) - num_commercialFEM => num_solver%get_asDict('Marc', defaultVal = emptyDict) + num_solver => config_numerics%get_dict('solver',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') diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index d6b0c7b77..9f326673b 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -125,7 +125,7 @@ subroutine FEM_utilities_init(num_mesh) 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 + write(petsc_optionsOrder,'(a,i0)') ' -mechFE_petscspace_degree ', p_s call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),err_PETSc) CHKERRQ(err_PETSc) From 99292c50819bf62b5e07bf85ce05f6209976378c Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 2 Aug 2023 16:02:51 +0200 Subject: [PATCH 04/10] prefix not required --- PRIVATE | 2 +- src/mesh/FEM_utilities.f90 | 15 +++++++-------- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/PRIVATE b/PRIVATE index 9a067cb27..12fe8685a 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 9a067cb2757df91b9242414779a8b52ac9e5ee02 +Subproject commit 12fe8685a2a2c1e37a8bb8e63e03130c99295023 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 9f326673b..cca975724 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -118,20 +118,19 @@ subroutine FEM_utilities_init(num_mesh) CHKERRQ(err_PETSc) CHKERRQ(err_PETSc) - 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' // & + petsc_options = misc_prefixOptions('-snes_type newtonls & + &-snes_linesearch_type cp -snes_ksp_ew & + &-snes_ksp_ew_rtol0 0.01 -snes_ksp_ew_rtolmax 0.01 & + &-ksp_type fgmres -ksp_max_it 25 ' // & num_mesh%get_asStr('PETSc_options',defaultVal=''), 'mechanical_') + + write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s + petsc_options = petsc_options // ' ' // petsc_optionsOrder 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) - CHKERRQ(err_PETSc) wgt = real(mesh_maxNips*mesh_NcpElemsGlobal,pREAL)**(-1) - end subroutine FEM_utilities_init From 8399a1170311c34558eadbd7a4e0b36d9d76c579 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 4 Sep 2023 11:39:45 -0400 Subject: [PATCH 05/10] more extensive numerics.yaml documentation --- examples/config/numerics.yaml | 110 +++++++++++++++++----------------- 1 file changed, 56 insertions(+), 54 deletions(-) diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index 487b6cf88..4ff1d3863 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -1,5 +1,5 @@ -# Available numerical parameters -# Case sensitive keys +# Default values of all available numerical parameters +# Please note that keys are case sensitive homogenization: mechanical: @@ -21,75 +21,77 @@ homogenization: solver: grid: - N_staggered_iter_max: 10 # max number of field level staggered iterations - 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 + N_cutback_max: 3 # maximum cut-back level (0: 1, 1: 0.5, 2: 0.25, etc) damage: - 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 + 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 + thermal: - N_iter_max: 100 # maximum iteration number - eps_abs_T: 1.0e-2 # absolute tolerance for thermal equilibrium - eps_rel_T: 1.0e-6 # relative tolerance for thermal equilibrium + N_iter_max: 100 # maximum iteration number + eps_abs_T: 1.0e-2 # absolute tolerance for thermal equilibrium + eps_rel_T: 1.0e-6 # relative tolerance for thermal equilibrium mechanical: - eps_abs_div(P): 1.0e-4 # absolute tolerance for fulfillment of stress equilibrium - eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium - eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC - 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) + N_iter_min: 1 # minimum iteration number + N_iter_max: 100 # maximum iteration number + eps_abs_div(P): 1.0e-4 # absolute tolerance for fulfillment of stress equilibrium + eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium + eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC + eps_rel_P: 1.0e-3 # relative tolerance for fulfillment of stress BC + update_gamma: false # update Gamma-operator with current dPdF (not possible if FFT:memory_efficient=true) FFT: - memory_efficient: True # Precalculate Gamma-operator (81 double per point) - divergence_correction: size+grid # Use size-independent divergence criterion - derivative: continuous # Approximation used for derivatives in Fourier space - FFTW_plan_mode: FFTW_MEASURE # planing-rigor flag, see manual on www.fftw.org - FFTW_timelimit: -1.0 # timelimit of plan creation for FFTW, see manual on www.fftw.org. -1.0: disable timelimit - PETSc_options: -snes_type ngmres -snes_ngmres_anderson # PETSc solver options - alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0. alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme - beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0. beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme - 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 + memory_efficient: true # precalculate Gamma-operator (81 doubles per point) + divergence_correction: size+grid # use size-independent divergence criterion {none, size, size+grid} + derivative: continuous # approximation used for derivatives in Fourier space {continuous, central_difference, FWBW_difference} + FFTW_plan_mode: FFTW_MEASURE # planning-rigor flags, see manual at https://www.fftw.org/fftw3_doc/Planner-Flags.html + FFTW_timelimit: -1.0 # time limit of plan creation for FFTW, see manual on www.fftw.org. (-1.0: disable time limit) + PETSc_options: -snes_type ngmres -snes_ngmres_anderson # PETSc solver options + alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0 (alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme) + beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0 ( beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme) + 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) + 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 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 + integrationorder: 2 # order of quadrature rule required + 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 phase: mechanical: - r_cutback_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation - r_cutback: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) - r_increase: 1.5 # factor to increase size of next step when previous step converged in phase state calculation - eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) - N_iter_state_max: 10 # state loop limit + r_cutback_min: 1.0e-3 # minimum (relative) size of step allowed during cutback in phase state calculation + r_cutback: 0.25 # factor to decrease size of step when cutback introduced in phase state calculation (value between 0 and 1) + r_increase: 1.5 # factor to increase size of next step when previous step converged in phase state calculation + eps_rel_state: 1.0e-6 # relative tolerance in phase state loop (abs tol provided by constitutive law) + N_iter_state_max: 10 # state loop limit plastic: - r_linesearch_Lp: 0.5 # factor to decrease the step due to non-convergence in Lp calculation - eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum - eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum - N_iter_Lp_max: 40 # stress loop limit for Lp - f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp - integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) + r_linesearch_Lp: 0.5 # factor to decrease the step due to non-convergence in Lp calculation + eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum + eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum + N_iter_Lp_max: 40 # stress loop limit for Lp + f_update_jacobi_Lp: 1 # frequency of Jacobian update of residuum in Lp + integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) + eigen: - r_linesearch_Li: 0.5 # factor to decrease the step due to non-convergence in Li calculation - eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum - eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum - N_iter_Li_max: 40 # stress loop limit for Li - f_update_jacobi_Li: 1 # frequency of Jacobian update of residuum in Li + r_linesearch_Li: 0.5 # factor to decrease the step due to non-convergence in Li calculation + eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum + eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum + N_iter_Li_max: 40 # stress loop limit for Li + f_update_jacobi_Li: 1 # frequency of updating the Jacobian of residuum in Li commercialFEM: - unitlength: 1 # physical length of one computational length unit + unitlength: 1 # physical length of one computational length unit generic: - random_seed: 0 # fixed seeding for pseudo-random number generator, Default 0: use random seed. - phi_min: 1.0e-6 # non-zero residual damage. + random_seed: 0 # fixed seeding for pseudo-random number generator (0: use random seed) + phi_min: 1.0e-6 # non-zero residual damage From 2434c322cadbb9a7a47fc28fc5f3269f509e5d3d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 5 Sep 2023 13:14:20 -0400 Subject: [PATCH 06/10] more sensible sorting of sections; additional polish --- examples/config/numerics.yaml | 68 +++++++++++++++++------------------ 1 file changed, 34 insertions(+), 34 deletions(-) diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index 4ff1d3863..e04888055 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -1,28 +1,10 @@ # Default values of all available numerical parameters # Please note that keys are case sensitive -homogenization: - mechanical: - RGC: - eps_abs_P: 1.0e+4 # absolute tolerance of RGC residuum (in Pa) - eps_rel_P: 1.0e-3 # relative ... - eps_abs_max: 1.0e+10 # absolute upper-limit of RGC residuum (in Pa) - eps_rel_max: 1.0e+2 # relative ... - Delta_a: 1.0e-7 # perturbation for computing penalty tangent - relevant_mismatch: 1.0e-5 # minimum threshold of mismatch - viscosity_exponent: 1.0e+0 # power (sensitivity rate) of numerical viscosity in RGC scheme - viscosity_modulus: 0.0e+0 # stress modulus of RGC numerical viscosity (zero = without numerical viscosity) - # suggestion: larger than the aTol_RGC but still far below the expected flow stress of material - dot_a_ref: 1.0e-3 # reference rate of relaxation (about the same magnitude as straining rate, possibly a bit higher) - dot_a_max: 1.0e+0 # threshold of maximum relaxation vector increment (if exceed this then cutback) - Delta_V_max: 1.0e-5 # maximum allowable relative volume discrepancy - Delta_V_modulus: 1.0e+12 - Delta_V_exponent: 5.0 - solver: grid: N_staggered_iter_max: 10 # max number of field-level staggered iterations - N_cutback_max: 3 # maximum cut-back level (0: 1, 1: 0.5, 2: 0.25, etc) + N_cutback_max: 3 # maximum cutback level (0: 1, 1: 0.5, 2: 0.25, etc) damage: N_iter_max: 100 # maximum iteration number @@ -41,7 +23,7 @@ solver: eps_rel_div(P): 5.0e-4 # relative tolerance for fulfillment of stress equilibrium eps_abs_P: 1.0e3 # absolute tolerance for fulfillment of stress BC eps_rel_P: 1.0e-3 # relative tolerance for fulfillment of stress BC - update_gamma: false # update Gamma-operator with current dPdF (not possible if FFT:memory_efficient=true) + update_gamma: false # update Gamma-operator with current dPdF (not possible if FFT: memory_efficient == true) FFT: memory_efficient: true # precalculate Gamma-operator (81 doubles per point) @@ -50,21 +32,28 @@ solver: FFTW_plan_mode: FFTW_MEASURE # planning-rigor flags, see manual at https://www.fftw.org/fftw3_doc/Planner-Flags.html FFTW_timelimit: -1.0 # time limit of plan creation for FFTW, see manual on www.fftw.org. (-1.0: disable time limit) PETSc_options: -snes_type ngmres -snes_ngmres_anderson # PETSc solver options - alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0 (alpha = 1.0 ==> AL scheme, alpha = 2.0 ==> accelerated scheme) - beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0 ( beta = 1.0 ==> AL scheme, beta = 2.0 ==> accelerated scheme) + alpha: 1.0 # polarization scheme parameter 0.0 < alpha < 2.0 (1.0: AL scheme, 2.0: accelerated scheme) + beta: 1.0 # polarization scheme parameter 0.0 < beta < 2.0 (1.0: AL scheme, 2.0: accelerated scheme) 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 - bbarstabilisation: false - integrationorder: 2 # order of quadrature rule required - 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 +homogenization: + mechanical: + RGC: + eps_abs_P: 1.0e+4 # absolute tolerance of RGC residuum (in Pa) + eps_rel_P: 1.0e-3 # relative ... + eps_abs_max: 1.0e+10 # absolute upper-limit of RGC residuum (in Pa) + eps_rel_max: 1.0e+2 # relative ... + Delta_a: 1.0e-7 # perturbation for computing penalty tangent + relevant_mismatch: 1.0e-5 # minimum threshold of mismatch + viscosity_exponent: 1.0e+0 # power (sensitivity rate) of numerical viscosity in RGC scheme + viscosity_modulus: 0.0e+0 # stress modulus of RGC numerical viscosity (0: without numerical viscosity) + # suggestion: larger than the aTol_RGC but still far below the expected flow stress of material + dot_a_ref: 1.0e-3 # reference rate of relaxation (about the same magnitude as straining rate, possibly a bit higher) + dot_a_max: 1.0e+0 # threshold of maximum relaxation vector increment (if exceed this then cutback) + Delta_V_max: 1.0e-5 # maximum allowable relative volume discrepancy + Delta_V_modulus: 1.0e+12 + Delta_V_exponent: 5.0 phase: mechanical: @@ -75,7 +64,7 @@ phase: N_iter_state_max: 10 # state loop limit plastic: - r_linesearch_Lp: 0.5 # factor to decrease the step due to non-convergence in Lp calculation + r_linesearch_Lp: 0.5 # factor to decrease the step if Lp calculation fails to converge eps_rel_Lp: 1.0e-6 # relative tolerance in Lp residuum eps_abs_Lp: 1.0e-8 # absolute tolerance in Lp residuum N_iter_Lp_max: 40 # stress loop limit for Lp @@ -83,7 +72,7 @@ phase: integrator_state: FPI # integration method (FPI = Fixed Point Iteration, Euler = Euler, AdaptiveEuler = Adaptive Euler, RK4 = classical 4th order Runge-Kutta, RKCK45 = 5th order Runge-Kutta Cash-Karp) eigen: - r_linesearch_Li: 0.5 # factor to decrease the step due to non-convergence in Li calculation + r_linesearch_Li: 0.5 # factor to decrease the step if Li calculation fails to converge eps_rel_Li: 1.0e-6 # relative tolerance in Li residuum eps_abs_Li: 1.0e-8 # absolute tolerance in Li residuum N_iter_Li_max: 40 # stress loop limit for Li @@ -95,3 +84,14 @@ commercialFEM: generic: random_seed: 0 # fixed seeding for pseudo-random number generator (0: use random seed) phi_min: 1.0e-6 # non-zero residual damage + +mesh: + maxCutBack: 3 # maximum cutback 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 + bbarstabilisation: false + integrationorder: 2 # order of quadrature rule required + 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 From e101b5cca8e18a805c65071b41718007bed2f9b2 Mon Sep 17 00:00:00 2001 From: Sharan Date: Wed, 6 Sep 2023 11:17:39 +0200 Subject: [PATCH 07/10] correct pointer --- src/mesh/discretization_mesh.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 8cecfbfb2..7790a82fd 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -101,7 +101,7 @@ subroutine discretization_mesh_init(restart) !-------------------------------------------------------------------------------- ! read numerics parameter num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict) - num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict) + num_mesh => num_solver%get_dict('mesh',defaultVal=emptyDict) p_i = num_mesh%get_asInt('p_i',defaultVal = 2) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16) From 252128a8638a39dcee601ee449516fef75dcbd20 Mon Sep 17 00:00:00 2001 From: Sharan Date: Wed, 6 Sep 2023 12:16:58 +0200 Subject: [PATCH 08/10] updated private repo --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index d85026fcc..1f8e4c450 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit d85026fccb1cad572a24ef19eefc36365df78a0f +Subproject commit 1f8e4c45073b0bd79f95caf034e993c72a74410d From 9c4fa16c7d8f117e6281306d36c2c30b2a9a6577 Mon Sep 17 00:00:00 2001 From: Sharan Date: Wed, 6 Sep 2023 19:57:40 +0200 Subject: [PATCH 09/10] new mech dict --- PRIVATE | 2 +- src/mesh/FEM_utilities.f90 | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/PRIVATE b/PRIVATE index 1f8e4c450..23a19f513 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 1f8e4c45073b0bd79f95caf034e993c72a74410d +Subproject commit 23a19f513e760e0b6d93ebb3bbf03ca3171c71d6 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 2389721b4..4480e412d 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -95,6 +95,9 @@ subroutine FEM_utilities_init(num_mesh) type(tDict), pointer, intent(in) :: & num_mesh + + type(tDict), pointer :: & + num_mech character(len=pSTRLEN) :: petsc_optionsOrder character(len=:), allocatable :: & petsc_options @@ -106,6 +109,8 @@ subroutine FEM_utilities_init(num_mesh) print'(/,1x,a)', '<<<+- FEM_utilities init -+>>>' + num_mech => num_mesh%get_dict('mechanical', defaultVal=emptyDict) + p_s = num_mesh%get_asInt('p_s',defaultVal = 2) p_i = num_mesh%get_asInt('p_i',defaultVal = p_s) @@ -123,7 +128,7 @@ subroutine FEM_utilities_init(num_mesh) &-snes_linesearch_type cp -snes_ksp_ew & &-snes_ksp_ew_rtol0 0.01 -snes_ksp_ew_rtolmax 0.01 & &-ksp_type fgmres -ksp_max_it 25 ' // & - num_mesh%get_asStr('PETSc_options',defaultVal=''), 'mechanical_') + num_mech%get_asStr('PETSc_options',defaultVal=''), 'mechanical_') write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s petsc_options = petsc_options // ' ' // petsc_optionsOrder From af955df891688949bb7062dce3a049e757538c7c Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 6 Sep 2023 15:11:48 -0400 Subject: [PATCH 10/10] added phi_min back to example; more consistency in code syntax --- examples/config/numerics.yaml | 3 ++- src/Marc/discretization_Marc.f90 | 4 ++-- src/mesh/DAMASK_mesh.f90 | 4 ++-- src/mesh/discretization_mesh.f90 | 2 +- src/mesh/mesh_mech_FEM.f90 | 16 ++++++++-------- 5 files changed, 15 insertions(+), 14 deletions(-) diff --git a/examples/config/numerics.yaml b/examples/config/numerics.yaml index 33bd267fc..cdc487589 100644 --- a/examples/config/numerics.yaml +++ b/examples/config/numerics.yaml @@ -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 diff --git a/src/Marc/discretization_Marc.f90 b/src/Marc/discretization_Marc.f90 index a5094b7fa..5e42180de 100644 --- a/src/Marc/discretization_Marc.f90 +++ b/src/Marc/discretization_Marc.f90 @@ -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) diff --git a/src/mesh/DAMASK_mesh.f90 b/src/mesh/DAMASK_mesh.f90 index 133a98997..370dc3535 100644 --- a/src/mesh/DAMASK_mesh.f90 +++ b/src/mesh/DAMASK_mesh.f90 @@ -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) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 7790a82fd..be9be3b19 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -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) diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index 403372db7..c566990e2 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -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), &