Merge branch 'restructure-numerics' into 'development'

Restructure numerics

Closes #271 and #63

See merge request damask/DAMASK!790

PRIVATE

@@ -1 +1 @@
-Subproject commit cea3296f78f5b48ab949364a80b738d3a7cf2194
+Subproject commit 23a19f513e760e0b6d93ebb3bbf03ca3171c71d6

examples/config/numerics.yaml

@@ -1,5 +1,58 @@
-# Available numerical parameters
+# Default values of all available numerical parameters
-# Case sensitive keys
+# Please note that keys are case sensitive
+
+solver:
+  grid:
+    N_staggered_iter_max: 10            # max number of field-level staggered iterations
+    N_cutback_max:        3             # maximum cutback 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
+      phi_min: 1.0e-6                   # residual integrity
+
+    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
+
+    mechanical:
+      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 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 (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:
+    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
+    p_i:                  2             # order of quadrature rule required
+    bbarstabilization:    false
+
+    mechanical:
+      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
+
 
 homogenization:
   mechanical:
@@ -11,7 +64,7 @@ homogenization:
       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)
+      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)
@@ -19,51 +72,6 @@ homogenization:
       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)
-
-    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
-    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
-
-    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)
-
-    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
-
-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
 
 phase:
   mechanical:
@@ -74,22 +82,19 @@ 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
       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
+      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
-      f_update_jacobi_Li:      1            # frequency of Jacobian update of residuum in 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
 
 generic:
-  random_seed:                 0            # fixed seeding for pseudo-random number generator, Default 0: use random seed.
+  random_seed:                 0        # fixed seeding for pseudo-random number generator (0: use random seed)
-  phi_min:                     1.0e-6       # non-zero residual damage.

src/Marc/discretization_Marc.f90

@@ -69,14 +69,16 @@ 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)
+  num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
-  mesh_unitlength = num_commercialFEM%get_asReal('unitlength',defaultVal=1.0_pREAL)                 ! set physical extent of a length unit in mesh
+  num_commercialFEM => num_solver%get_dict('Marc',defaultVal=emptyDict)
-  if (mesh_unitlength <= 0.0_pREAL) call IO_error(301,'unitlength')
+  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,'unit_length')
 
   call inputRead(elem,node0_elem,connectivity_elem,materialAt)
   nElems = size(connectivity_elem,2)

src/grid/grid_damage_spectral.f90

@@ -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))

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)
+  num_solver => config_numerics%get_dict('solver',defaultVal=emptyDict)
-  stagItMax  = num_mesh%get_asInt('maxStaggeredIter',defaultVal=10)
+  num_mesh   => num_solver%get_dict('mesh',defaultVal=emptyDict)
-  maxCutBack = num_mesh%get_asInt('maxCutBack',defaultVal=3)
+  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 (stagItMax < 0)  call IO_error(301,ext_msg='N_staggered_iter_max')
-  if (maxCutBack < 0) call IO_error(301,ext_msg='maxCutBack')
+  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_Utilities_init(num_mesh)
-  call FEM_mechanical_init(loadCases(1)%fieldBC(1))
+  call FEM_mechanical_init(loadCases(1)%fieldBC(1),num_mesh)
   call config_numerics_deallocate()
 
   if (worldrank == 0) then

src/mesh/FEM_utilities.f90

@@ -16,6 +16,7 @@ module FEM_utilities
   use prec
   use config
   use math
+  use misc
   use IO
   use discretization_mesh
   use homogenization
@@ -90,11 +91,16 @@ contains
 !--------------------------------------------------------------------------------------------------
 !> @brief Allocate all neccessary fields.
 !--------------------------------------------------------------------------------------------------
-subroutine FEM_utilities_init
+subroutine FEM_utilities_init(num_mesh)
 
-  character(len=pSTRLEN) :: petsc_optionsOrder
+  type(tDict), pointer, intent(in) :: &
-  type(tDict), pointer :: &
     num_mesh
+
+  type(tDict), pointer :: &
+    num_mech
+  character(len=pSTRLEN) :: petsc_optionsOrder
+  character(len=:), allocatable :: &
+    petsc_options
   integer :: &
     p_s, &                                                                                          !< order of shape functions
     p_i                                                                                             !< integration order (quadrature rule)
@@ -103,7 +109,7 @@ subroutine FEM_utilities_init
 
   print'(/,1x,a)',   '<<<+-  FEM_utilities init  -+>>>'
 
-  num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict)
+  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)
@@ -117,20 +123,20 @@ 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 &
+
-                               &-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew &
+  petsc_options = misc_prefixOptions('-snes_type newtonls &
-                               &-mechanical_snes_ksp_ew_rtol0 0.01 -mechanical_snes_ksp_ew_rtolmax 0.01 &
+                               &-snes_linesearch_type cp -snes_ksp_ew &
-                               &-mechanical_ksp_type fgmres -mechanical_ksp_max_it 25', err_PETSc)
+                               &-snes_ksp_ew_rtol0 0.01 -snes_ksp_ew_rtolmax 0.01 &
-  CHKERRQ(err_PETSc)
+                               &-ksp_type fgmres -ksp_max_it 25 ' // &
-  call PetscOptionsInsertString(PETSC_NULL_OPTIONS,num_mesh%get_asStr('PETSc_options',defaultVal=''),err_PETSc)
+                                num_mech%get_asStr('PETSc_options',defaultVal=''), 'mechanical_')
-  CHKERRQ(err_PETSc)
+
   write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', p_s
-  call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),err_PETSc)
+  petsc_options = petsc_options // ' ' // petsc_optionsOrder
+  call PetscOptionsInsertString(PETSC_NULL_OPTIONS,petsc_options,err_PETSc)
   CHKERRQ(err_PETSc)
 
   wgt = real(mesh_maxNips*mesh_NcpElemsGlobal,pREAL)**(-1)
 
-
 end subroutine FEM_utilities_init
 
 

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   => 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)

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
@@ -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(tDict), pointer, intent(in)       :: num_mesh
 
   DM                                     :: mechanical_mesh
   PetscFE                                :: mechFE
@@ -126,23 +127,24 @@ subroutine FEM_mechanical_init(fieldBC)
   character(len=*), parameter            :: prefix = 'mechFE_'
   PetscErrorCode                         :: err_PETSc
   real(pREAL), dimension(3,3) :: devNull
-  type(tDict), pointer :: &
+  type(tDict), pointer                   :: num_mech
-    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_mech => num_mesh%get_dict('mechanical', 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%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)
 
-  if (num%itmax <= 1)                       call IO_error(301,ext_msg='itmax')
+  num%p_i               = int(num_mesh%get_asInt('p_i',defaultVal=2),pPETSCINT)
-  if (num%eps_struct_rtol <= 0.0_pREAL)     call IO_error(301,ext_msg='eps_struct_rtol')
+  num%BBarStabilization = num_mesh%get_asBool('bbarstabilization',defaultVal=.false.)
-  if (num%eps_struct_atol <= 0.0_pREAL)     call IO_error(301,ext_msg='eps_struct_atol')
+
+  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)
+
+  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
@@ -437,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
@@ -588,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))
@@ -604,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), &