Merge branch '260-remove-debug-yaml-assocciated-code-tests' into 'development'

debug.yaml causes more work than it saves

Closes #260

See merge request damask/DAMASK!717
This commit is contained in:
Franz Roters 2023-02-14 08:12:59 +00:00
commit b4e706f099
17 changed files with 29 additions and 288 deletions

View File

@ -1,12 +0,0 @@
phase: [basic, extensive, selective]
materialpoint: [basic, extensive, selective]
# options for selective debugging
element: 1
integrationpoint: 1
constituent: 1
# solver-specific
mesh: [PETSc]
grid: [basic, rotation, PETSc]
Marc: [basic]

View File

@ -130,7 +130,7 @@ subroutine CLI_init
print'(a)', ' Make sure the file "material.yaml" exists in the working' print'(a)', ' Make sure the file "material.yaml" exists in the working'
print'(a)', ' directory.' print'(a)', ' directory.'
print'(a)', ' For further configuration place "numerics.yaml"' print'(a)', ' For further configuration place "numerics.yaml"'
print'(a)',' and "debug.yaml" in that directory.' print'(a)',' in that directory.'
print'(/,a)',' --restart N' print'(/,a)',' --restart N'
print'(a)', ' Reads in increment N and continues with calculating' print'(a)', ' Reads in increment N and continues with calculating'
print'(a)', ' increment N+1, N+2, ... based on this.' print'(a)', ' increment N+1, N+2, ... based on this.'

View File

@ -484,8 +484,6 @@ subroutine IO_error(error_ID,ext_msg,label1,ID1,label2,ID2)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! user errors ! user errors
case (602)
msg = 'invalid selection for debug'
case (603) case (603)
msg = 'invalid data for table' msg = 'invalid data for table'

View File

@ -282,34 +282,15 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
logical, save :: & logical, save :: &
lastIncConverged = .false., & !< needs description lastIncConverged = .false., & !< needs description
outdatedByNewInc = .false., & !< needs description outdatedByNewInc = .false., & !< needs description
materialpoint_init_done = .false., & !< remember whether init has been done already materialpoint_init_done = .false. !< remember whether init has been done already
debug_basic = .true.
type(tList), pointer :: &
debug_Marc ! pointer to Marc debug options
if (debug_basic) then
print'(a,/,i8,i8,i2)', ' MSC.Marc information on shape of element(2), IP:', m, nn
print'(a,2(i1))', ' Jacobian: ', ngens,ngens
print'(a,i1)', ' Direct stress: ', ndi
print'(a,i1)', ' Shear stress: ', nshear
print'(a,i2)', ' DoF: ', ndeg
print'(a,i2)', ' Coordinates: ', ncrd
print'(a,i12)', ' Nodes: ', nnode
print'(a,i1)', ' Deformation gradient: ', itel
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
transpose(ffn)
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
transpose(ffn1)
end if
defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
call omp_set_num_threads(1_pI32) ! no openMP call omp_set_num_threads(1_pI32) ! no openMP
if (.not. materialpoint_init_done) then if (.not. materialpoint_init_done) then
materialpoint_init_done = .true. materialpoint_init_done = .true.
call materialpoint_initAll call materialpoint_initAll()
debug_Marc => config_debug%get_list('Marc',defaultVal=emptyList)
debug_basic = debug_Marc%contains('basic')
end if end if
computationMode = 0 ! save initialization value, since it does not result in any calculation computationMode = 0 ! save initialization value, since it does not result in any calculation

View File

@ -59,8 +59,7 @@ subroutine discretization_Marc_init
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
materialAt materialAt
integer:: & integer:: &
Nelems, & !< total number of elements in the mesh Nelems !< total number of elements in the mesh
debug_e, debug_i
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
IP_reshaped IP_reshaped
@ -75,9 +74,6 @@ subroutine discretization_Marc_init
print'(/,a)', ' <<<+- discretization_Marc init -+>>>'; flush(6) print'(/,a)', ' <<<+- discretization_Marc init -+>>>'; flush(6)
debug_e = config_debug%get_asInt('element',defaultVal=1)
debug_i = config_debug%get_asInt('integrationpoint',defaultVal=1)
num_commercialFEM => config_numerics%get_dict('commercialFEM',defaultVal = emptyDict) num_commercialFEM => config_numerics%get_dict('commercialFEM',defaultVal = emptyDict)
mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',defaultVal=1.0_pReal) ! set physical extent of a length unit in mesh mesh_unitlength = num_commercialFEM%get_asFloat('unitlength',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,'unitlength')
@ -85,9 +81,6 @@ subroutine discretization_Marc_init
call inputRead(elem,node0_elem,connectivity_elem,materialAt) call inputRead(elem,node0_elem,connectivity_elem,materialAt)
nElems = size(connectivity_elem,2) nElems = size(connectivity_elem,2)
if (debug_e < 1 .or. debug_e > nElems) call IO_error(602,'element')
if (debug_i < 1 .or. debug_i > elem%nIPs) call IO_error(602,'IP')
allocate(cellNodeDefinition(elem%nNodes-1)) allocate(cellNodeDefinition(elem%nNodes-1))
allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems)) allocate(connectivity_cell(elem%NcellNodesPerCell,elem%nIPs,nElems))

View File

@ -50,18 +50,6 @@ module materialpoint_Marc
type(tNumerics), private :: num type(tNumerics), private :: num
type, private :: tDebugOptions
logical :: &
basic, &
extensive, &
selective
integer:: &
element, &
ip
end type tDebugOptions
type(tDebugOptions), private :: debugmaterialpoint
public :: & public :: &
materialpoint_general, & materialpoint_general, &
materialpoint_initAll, & materialpoint_initAll, &
@ -99,36 +87,16 @@ end subroutine materialpoint_initAll
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocate the arrays defined in module materialpoint and initialize them !> @brief Allocate the arrays defined in module materialpoint and initialize them.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine materialpoint_init() subroutine materialpoint_init()
type(tList), pointer :: &
debug_materialpoint
print'(/,1x,a)', '<<<+- materialpoint init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- materialpoint init -+>>>'; flush(IO_STDOUT)
allocate(materialpoint_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal) allocate(materialpoint_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(materialpoint_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal) allocate(materialpoint_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(materialpoint_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal) allocate(materialpoint_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
!------------------------------------------------------------------------------
! read debug options
debug_materialpoint => config_debug%get_list('materialpoint',defaultVal=emptyList)
debugmaterialpoint%basic = debug_materialpoint%contains('basic')
debugmaterialpoint%extensive = debug_materialpoint%contains('extensive')
debugmaterialpoint%selective = debug_materialpoint%contains('selective')
debugmaterialpoint%element = config_debug%get_asInt('element',defaultVal = 1)
debugmaterialpoint%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
if (debugmaterialpoint%basic) then
print'(a32,1x,6(i8,1x))', 'materialpoint_cs: ', shape(materialpoint_cs)
print'(a32,1x,6(i8,1x))', 'materialpoint_dcsdE: ', shape(materialpoint_dcsdE)
print'(a32,1x,6(i8,1x),/)', 'materialpoint_dcsdE_knownGood: ', shape(materialpoint_dcsdE_knownGood)
flush(IO_STDOUT)
end if
end subroutine materialpoint_init end subroutine materialpoint_init
@ -164,16 +132,6 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
elCP = discretization_Marc_FEM2DAMASK_elem(elFE) elCP = discretization_Marc_FEM2DAMASK_elem(elFE)
ce = discretization_Marc_FEM2DAMASK_cell(ip,elFE) ce = discretization_Marc_FEM2DAMASK_cell(ip,elFE)
if (debugmaterialpoint%basic .and. elCP == debugmaterialpoint%element .and. ip == debugmaterialpoint%ip) then
print'(/,a)', '#############################################'
print'(a1,a22,1x,i8,a13)', '#','element', elCP, '#'
print'(a1,a22,1x,i8,a13)', '#','ip', ip, '#'
print'(a1,a22,1x,i8,a13)', '#','cycleCounter', cycleCounter, '#'
print'(a1,a22,1x,i8,a13)', '#','computationMode',mode, '#'
if (terminallyIll) &
print'(a,/)', '# --- terminallyIll --- #'
print'(a,/)', '#############################################'; flush (6)
end if
if (iand(mode, materialpoint_BACKUPJACOBIAN) /= 0) & if (iand(mode, materialpoint_BACKUPJACOBIAN) /= 0) &
materialpoint_dcsde_knownGood = materialpoint_dcsde materialpoint_dcsde_knownGood = materialpoint_dcsde
@ -194,7 +152,6 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6) materialpoint_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
else validCalculation else validCalculation
if (debugmaterialpoint%extensive) print'(a,i8,1x,i2)', '<< materialpoint >> calculation for elFE ip ',elFE,ip
call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip) call homogenization_mechanical_response(dt,(elCP-1)*discretization_nIPs + ip,(elCP-1)*discretization_nIPs + ip)
if (.not. terminallyIll) & if (.not. terminallyIll) &
call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP]) call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
@ -232,15 +189,6 @@ subroutine materialpoint_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip,
end if terminalIllness end if terminalIllness
end if validCalculation end if validCalculation
if (debugmaterialpoint%extensive &
.and. ((debugmaterialpoint%element == elCP .and. debugmaterialpoint%ip == ip) .or. .not. debugmaterialpoint%selective)) then
print'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)', &
'<< materialpoint >> stress/MPa at elFE ip ', elFE, ip, materialpoint_cs(1:6,ip,elCP)*1.0e-6_pReal
print'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))', &
'<< materialpoint >> Jacobian/GPa at elFE ip ', elFE, ip, transpose(materialpoint_dcsdE(1:6,1:6,ip,elCP))*1.0e-9_pReal
flush(IO_STDOUT)
end if
end if end if
if (all(abs(materialpoint_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) & if (all(abs(materialpoint_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) &

View File

@ -1,6 +1,6 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Read in the configuration of material, numerics, and debug from their respective file !> @brief Read in the material and numerics configuration from their respective file.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module config module config
use IO use IO
@ -14,8 +14,7 @@ module config
type(tDict), pointer, public :: & type(tDict), pointer, public :: &
config_material, & config_material, &
config_numerics, & config_numerics
config_debug
public :: & public :: &
config_init, & config_init, &
@ -24,7 +23,7 @@ module config
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Real *.yaml configuration files. !> @brief Read *.yaml configuration files.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine config_init() subroutine config_init()
@ -32,7 +31,6 @@ subroutine config_init()
call parse_material() call parse_material()
call parse_numerics() call parse_numerics()
call parse_debug()
end subroutine config_init end subroutine config_init
@ -95,41 +93,9 @@ subroutine parse_numerics()
end subroutine parse_numerics end subroutine parse_numerics
!--------------------------------------------------------------------------------------------------
!> @brief Read debug.yaml.
!--------------------------------------------------------------------------------------------------
subroutine parse_debug()
logical :: fileExists
character(len=:), allocatable :: fileContent
config_debug => emptyDict
inquire(file='debug.yaml', exist=fileExists)
if (fileExists) then
if (worldrank == 0) then
print'(1x,a)', 'reading debug.yaml'; flush(IO_STDOUT)
fileContent = IO_read('debug.yaml')
if (len(fileContent) > 0) then
call result_openJobFile(parallel=.false.)
call result_writeDataset_str(fileContent,'setup','debug.yaml','debug configuration')
call result_closeJobFile
end if
end if
call parallelization_bcast_str(fileContent)
config_debug => YAML_parse_str_asDict(fileContent)
end if
end subroutine parse_debug
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Deallocate config_material. !> @brief Deallocate config_material.
!ToDo: deallocation of numerics and debug (optional) !ToDo: deallocation of numerics (optional)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine config_deallocate subroutine config_deallocate

View File

@ -65,8 +65,7 @@ subroutine discretization_grid_init(restart)
materialAt, materialAt_global materialAt, materialAt_global
integer :: & integer :: &
j, & j
debug_element, debug_ip
integer(MPI_INTEGER_KIND) :: err_MPI integer(MPI_INTEGER_KIND) :: err_MPI
integer(C_INTPTR_T) :: & integer(C_INTPTR_T) :: &
devNull, z, z_offset devNull, z, z_offset
@ -163,13 +162,6 @@ subroutine discretization_grid_init(restart)
call geometry_plastic_nonlocal_setIPareaNormal (cellSurfaceNormal(product(myGrid))) call geometry_plastic_nonlocal_setIPareaNormal (cellSurfaceNormal(product(myGrid)))
call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(myGrid)) call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(myGrid))
!-------------------------------------------------------------------------------------------------
! debug parameters
debug_element = config_debug%get_asInt('element',defaultVal=1)
if (debug_element < 1 .or. debug_element > product(myGrid)) call IO_error(602,ext_msg='element')
debug_ip = config_debug%get_asInt('integrationpoint',defaultVal=1)
if (debug_ip /= 1) call IO_error(602,ext_msg='IP')
end subroutine discretization_grid_init end subroutine discretization_grid_init

View File

@ -50,8 +50,6 @@ module grid_mechanical_FEM
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
logical :: debugRotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: mechanical_grid DM :: mechanical_grid
@ -121,19 +119,12 @@ subroutine grid_mechanical_FEM_init
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid num_grid
type(tList), pointer :: &
debug_grid
character(len=pStringLen) :: & character(len=pStringLen) :: &
extmsg = '' extmsg = ''
print'(/,1x,a)', '<<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- grid_mechanical_FEM init -+>>>'; flush(IO_STDOUT)
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => config_debug%get_list('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict)
@ -565,7 +556,8 @@ subroutine formResidual(da_local,x_local, &
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter+1, '≤', num%itmax
if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & if (any(dNeq(params%rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pReal)))) &
print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
'deformation gradient aim =', transpose(F_aim) 'deformation gradient aim =', transpose(F_aim)

View File

@ -49,8 +49,6 @@ module grid_mechanical_spectral_basic
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
logical :: debugRotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: da DM :: da
@ -117,8 +115,6 @@ subroutine grid_mechanical_spectral_basic_init()
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid num_grid
type(tList), pointer :: &
debug_grid
character(len=pStringLen) :: & character(len=pStringLen) :: &
extmsg = '' extmsg = ''
@ -131,11 +127,6 @@ subroutine grid_mechanical_spectral_basic_init()
print'( 1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:3145, 2015' print'( 1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:3145, 2015'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006' print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => config_debug%get_list('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict)
@ -518,7 +509,8 @@ subroutine formResidual(residual_subdomain, F, &
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & if (any(dNeq(params%rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pReal)))) &
print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
'deformation gradient aim =', transpose(F_aim) 'deformation gradient aim =', transpose(F_aim)

View File

@ -54,8 +54,6 @@ module grid_mechanical_spectral_polarisation
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
logical :: debugRotation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! PETSc data ! PETSc data
DM :: da DM :: da
@ -130,8 +128,6 @@ subroutine grid_mechanical_spectral_polarisation_init()
integer(HID_T) :: fileHandle, groupHandle integer(HID_T) :: fileHandle, groupHandle
type(tDict), pointer :: & type(tDict), pointer :: &
num_grid num_grid
type(tList), pointer :: &
debug_grid
character(len=pStringLen) :: & character(len=pStringLen) :: &
extmsg = '' extmsg = ''
@ -141,10 +137,6 @@ subroutine grid_mechanical_spectral_polarisation_init()
print'(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:3145, 2015' print'(/,1x,a)', 'P. Shanthraj et al., International Journal of Plasticity 66:3145, 2015'
print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006' print'( 1x,a)', 'https://doi.org/10.1016/j.ijplas.2014.02.006'
!-------------------------------------------------------------------------------------------------
! debugging options
debug_grid => config_debug%get_list('grid',defaultVal=emptyList)
debugRotation = debug_grid%contains('rotation')
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! read numerical parameters and do sanity checks ! read numerical parameters and do sanity checks
@ -595,7 +587,8 @@ subroutine formResidual(residual_subdomain, FandF_tau, &
newIteration: if (totalIter <= PETScIter) then newIteration: if (totalIter <= PETScIter) then
totalIter = totalIter + 1 totalIter = totalIter + 1
print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax print'(1x,a,3(a,i0))', trim(incInfo), ' @ Iteration ', num%itmin, '≤',totalIter, '≤', num%itmax
if (debugRotation) print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & if (any(dNeq(params%rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pReal)))) &
print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.)) 'deformation gradient aim (lab) =', transpose(params%rotation_BC%rotate(F_aim,active=.true.))
print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', & print'(/,1x,a,/,2(3(f12.7,1x)/),3(f12.7,1x))', &
'deformation gradient aim =', transpose(F_aim) 'deformation gradient aim =', transpose(F_aim)

View File

@ -64,13 +64,6 @@ module spectral_utilities
planScalarForth, & !< FFTW MPI plan s(x) to s(k) planScalarForth, & !< FFTW MPI plan s(x) to s(k)
planScalarBack !< FFTW MPI plan s(k) to s(x) planScalarBack !< FFTW MPI plan s(k) to s(x)
!--------------------------------------------------------------------------------------------------
! variables controlling debugging
logical :: &
debugGeneral, & !< general debugging of spectral solver
debugRotation, & !< also printing out results in lab frame
debugPETSc !< use some in debug defined options for more verbose PETSc solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type, public :: tSolutionState !< return type of solution from spectral solver variants type, public :: tSolutionState !< return type of solution from spectral solver variants
@ -131,12 +124,7 @@ module spectral_utilities
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags, create plans for FFTW !> @brief Allocate all neccessary fields and create plans for FFTW.
!> @details Sets the debug levels for general, divergence, restart, and FFTW from the bitwise coding
!> provided by the debug module to logicals.
!> Allocate all fields used by FFTW and create the corresponding plans depending on the debug
!> level chosen.
!> Initializes FFTW.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine spectral_utilities_init() subroutine spectral_utilities_init()
@ -157,12 +145,8 @@ subroutine spectral_utilities_init()
integer(C_INTPTR_T), parameter :: & integer(C_INTPTR_T), parameter :: &
vectorSize = 3_C_INTPTR_T, & vectorSize = 3_C_INTPTR_T, &
tensorSize = 9_C_INTPTR_T tensorSize = 9_C_INTPTR_T
character(len=*), parameter :: &
PETSCDEBUG = ' -snes_view -snes_monitor '
type(tDict) , pointer :: & type(tDict) , pointer :: &
num_grid num_grid
type(tList) , pointer :: &
debug_grid
print'(/,1x,a)', '<<<+- spectral_utilities init -+>>>' print'(/,1x,a)', '<<<+- spectral_utilities init -+>>>'
@ -179,25 +163,10 @@ subroutine spectral_utilities_init()
print'( 1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019' print'( 1x,a)', 'P. Shanthraj et al., Handbook of Mechanics of Materials, 2019'
print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80' print'( 1x,a)', 'https://doi.org/10.1007/978-981-10-6855-3_80'
!--------------------------------------------------------------------------------------------------
! set debugging parameters
num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict) num_grid => config_numerics%get_dict('grid',defaultVal=emptyDict)
debug_grid => config_debug%get_List('grid',defaultVal=emptyList)
debugGeneral = debug_grid%contains('basic')
debugRotation = debug_grid%contains('rotation')
debugPETSc = debug_grid%contains('PETSc')
if (debugPETSc) print'(3(/,1x,a),/)', &
'Initializing PETSc with debug options: ', &
trim(PETScDebug), &
'add more using the "PETSc_options" keyword in numerics.yaml'
flush(IO_STDOUT)
call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc) call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),err_PETSc)
CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,& call PetscOptionsInsertString(PETSC_NULL_OPTIONS,&
num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc) num_grid%get_asString('PETSc_options',defaultVal=''),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -704,13 +673,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
if (size_reduced > 0) then if (size_reduced > 0) then
temp99_real = math_3333to99(rot_BC%rotate(C)) temp99_real = math_3333to99(rot_BC%rotate(C))
if (debugGeneral) then
print'(/,1x,a)', '... updating masked compliance ............................................'
print'(/,1x,a,/,8(9(2x,f12.7,1x)/),9(2x,f12.7,1x))', &
'Stiffness C (load) / GPa =', transpose(temp99_Real)*1.0e-9_pReal
flush(IO_STDOUT)
end if
do i = 1,9; do j = 1,9 do i = 1,9; do j = 1,9
mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j) mask(i,j) = mask_stressVector(i) .and. mask_stressVector(j)
end do; end do end do; end do
@ -724,7 +686,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
! check if inversion was successful ! check if inversion was successful
sTimesC = matmul(c_reduced,s_reduced) sTimesC = matmul(c_reduced,s_reduced)
errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal)) errmatinv = errmatinv .or. any(dNeq(sTimesC,math_eye(size_reduced),1.0e-12_pReal))
if (debugGeneral .or. errmatinv) then if (errmatinv) then
write(formatString, '(i2)') size_reduced write(formatString, '(i2)') size_reduced
formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))' formatString = '(/,1x,a,/,'//trim(formatString)//'('//trim(formatString)//'(2x,es9.2,1x)/))'
print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced)) print trim(formatString), 'C * S (load) ', transpose(matmul(c_reduced,s_reduced))
@ -738,12 +700,6 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
utilities_maskedCompliance = math_99to3333(temp99_Real) utilities_maskedCompliance = math_99to3333(temp99_Real)
if (debugGeneral) then
print'(/,1x,a,/,9(9(2x,f10.5,1x)/),9(2x,f10.5,1x))', &
'Masked Compliance (load) * GPa =', transpose(temp99_Real)*1.0e9_pReal
flush(IO_STDOUT)
end if
end function utilities_maskedCompliance end function utilities_maskedCompliance
@ -825,9 +781,12 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,P_av,9_MPI_INTEGER_KIND,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,err_MPI)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (debugRotation) print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & if (present(rotation_BC)) then
if (any(dNeq(rotation_BC%asQuaternion(), real([1.0, 0.0, 0.0, 0.0],pReal)))) &
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal 'Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal
if (present(rotation_BC)) P_av = rotation_BC%rotate(P_av) P_av = rotation_BC%rotate(P_av)
end if
print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', & print'(/,1x,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
'Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pReal 'Piola--Kirchhoff stress / MPa =', transpose(P_av)*1.e-6_pReal
flush(IO_STDOUT) flush(IO_STDOUT)

View File

@ -88,21 +88,17 @@ contains
!ToDo: use functions in variable call !ToDo: use functions in variable call
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, sets debug flags !> @brief Allocate all neccessary fields.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_utilities_init subroutine FEM_utilities_init
character(len=pStringLen) :: petsc_optionsOrder character(len=pStringLen) :: petsc_optionsOrder
type(tDict), pointer :: & type(tDict), pointer :: &
num_mesh, & num_mesh
debug_mesh ! pointer to mesh debug options
integer :: & integer :: &
p_s, & !< order of shape functions p_s, & !< order of shape functions
p_i !< integration order (quadrature rule) p_i !< integration order (quadrature rule)
character(len=*), parameter :: &
PETSCDEBUG = ' -snes_view -snes_monitor '
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
logical :: debugPETSc !< use some in debug defined options for more verbose PETSc solution
print'(/,1x,a)', '<<<+- FEM_utilities init -+>>>' print'(/,1x,a)', '<<<+- FEM_utilities init -+>>>'
@ -117,17 +113,9 @@ subroutine FEM_utilities_init
if (p_i < max(1,p_s-1) .or. p_i > p_s) & if (p_i < max(1,p_s-1) .or. p_i > p_s) &
call IO_error(821,ext_msg='integration order (p_i) out of bounds') call IO_error(821,ext_msg='integration order (p_i) out of bounds')
debug_mesh => config_debug%get_dict('mesh',defaultVal=emptyDict)
debugPETSc = debug_mesh%contains('PETSc')
if (debugPETSc) print'(3(/,1x,a),/)', &
'Initializing PETSc with debug options: ', &
trim(PETScDebug), &
'add more using the "PETSc_options" keyword in numerics.yaml'
flush(IO_STDOUT) flush(IO_STDOUT)
call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc) call PetscOptionsClear(PETSC_NULL_OPTIONS,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
if (debugPETSc) call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(PETSCDEBUG),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls & call PetscOptionsInsertString(PETSC_NULL_OPTIONS,'-mechanical_snes_type newtonls &
&-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew & &-mechanical_snes_linesearch_type cp -mechanical_snes_ksp_ew &

View File

@ -78,8 +78,7 @@ subroutine discretization_mesh_init(restart)
PetscInt :: dimPlex, & PetscInt :: dimPlex, &
mesh_Nnodes, & !< total number of nodes in mesh mesh_Nnodes, & !< total number of nodes in mesh
j, & j
debug_element, debug_ip
PetscSF :: sf PetscSF :: sf
DM :: globalMesh DM :: globalMesh
PetscInt :: nFaceSets, Nboundaries, NelemsGlobal, Nelems PetscInt :: nFaceSets, Nboundaries, NelemsGlobal, Nelems
@ -103,11 +102,6 @@ subroutine discretization_mesh_init(restart)
num_mesh => config_numerics%get_dict('mesh',defaultVal=emptyDict) num_mesh => config_numerics%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)
!---------------------------------------------------------------------------------
! read debug parameters
debug_element = config_debug%get_asInt('element',defaultVal=1)
debug_ip = config_debug%get_asInt('integrationpoint',defaultVal=1)
#if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16) #if (PETSC_VERSION_MAJOR==3 && PETSC_VERSION_MINOR>16)
call DMPlexCreateFromFile(PETSC_COMM_WORLD,CLI_geomFile,'n/a',PETSC_TRUE,globalMesh,err_PETSc) call DMPlexCreateFromFile(PETSC_COMM_WORLD,CLI_geomFile,'n/a',PETSC_TRUE,globalMesh,err_PETSc)
#else #else
@ -182,9 +176,6 @@ subroutine discretization_mesh_init(restart)
end do end do
materialAt = materialAt + 1_pPETSCINT materialAt = materialAt + 1_pPETSCINT
if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element')
if (debug_ip < 1 .or. debug_ip > mesh_maxNips) call IO_error(602,ext_msg='IP')
allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
mesh_node0(1:dimPlex,:) = reshape(mesh_node0_temp,[dimPlex,mesh_Nnodes]) mesh_node0(1:dimPlex,:) = reshape(mesh_node0_temp,[dimPlex,mesh_Nnodes])

View File

@ -76,17 +76,6 @@ module phase
type(tNumerics) :: num ! numerics parameters. Better name? type(tNumerics) :: num ! numerics parameters. Better name?
type :: tDebugOptions
logical :: &
basic, &
extensive, &
selective
integer :: &
element, &
ip, &
grain
end type tDebugOptions
type(tPlasticState), allocatable, dimension(:), public :: & type(tPlasticState), allocatable, dimension(:), public :: &
plasticState plasticState
type(tState), allocatable, dimension(:), public :: & type(tState), allocatable, dimension(:), public :: &
@ -334,7 +323,6 @@ module phase
end interface end interface
type(tDebugOptions) :: debugConstitutive
#if __INTEL_COMPILER >= 1900 #if __INTEL_COMPILER >= 1900
public :: & public :: &
prec, & prec, &
@ -390,21 +378,10 @@ subroutine phase_init
type(tDict), pointer :: & type(tDict), pointer :: &
phases, & phases, &
phase phase
type(tList), pointer :: &
debug_constitutive
print'(/,1x,a)', '<<<+- phase init -+>>>'; flush(IO_STDOUT) print'(/,1x,a)', '<<<+- phase init -+>>>'; flush(IO_STDOUT)
debug_constitutive => config_debug%get_list('phase', defaultVal=emptyList)
debugConstitutive%basic = debug_constitutive%contains('basic')
debugConstitutive%extensive = debug_constitutive%contains('extensive')
debugConstitutive%selective = debug_constitutive%contains('selective')
debugConstitutive%element = config_debug%get_asInt('element', defaultVal = 1)
debugConstitutive%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1)
debugConstitutive%grain = config_debug%get_asInt('constituent', defaultVal = 1)
phases => config_material%get_dict('phase') phases => config_material%get_dict('phase')
allocate(phase_lattice(phases%length)) allocate(phase_lattice(phases%length))
allocate(phase_cOverA(phases%length),source=-1.0_pReal) allocate(phase_cOverA(phases%length),source=-1.0_pReal)

View File

@ -1114,12 +1114,6 @@ module subroutine nonlocal_dotState(Mp,timestep, &
if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) & if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) &
.or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> evolution rate leads to negative density at ph ',ph,' en ',en
print'(a)', '<< CONST >> enforcing cutback !!!'
end if
#endif
plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
else else
dot%rho(:,en) = pack(rhoDot,.true.) dot%rho(:,en) = pack(rhoDot,.true.)
@ -1219,17 +1213,6 @@ function rhoDotFlux(timestep,ph,en)
if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ... if (any( abs(dot_gamma) > 0.0_pReal & ! any active slip system ...
.and. prm%C_CFL * abs(v0) * timestep & .and. prm%C_CFL * abs(v0) * timestep &
> geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) > geom(ph)%V_0(en)/ maxval(geom(ph)%IParea(:,en)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
#ifdef DEBUG
if (debugConstitutive%extensive) then
print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at ph ',ph,' en ',en
print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at ', &
maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
.and. prm%C_CFL * abs(v0) * timestep &
> geom(ph)%V_0(en) / maxval(geom(ph)%IParea(:,en))), &
' at a timestep of ',timestep
print*, '<< CONST >> enforcing cutback !!!'
end if
#endif
rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback rhoDotFlux = IEEE_value(1.0_pReal,IEEE_quiet_NaN) ! enforce cutback
return return
end if end if