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

Merge branch 'mesh-load-YAML' into 'development' 

Mesh load yaml

See merge request damask/DAMASK!876
This commit is contained in:
Martin Diehl 2023-12-21 18:38:36 +00:00
parent 1933c896b2 0518cc6730
commit af5bbed003
10 changed files with 213 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 29ef436acca5417aebc945b688642c34697af911 Subproject commit 62df7f24f2a95fda255f7d20b130afcfeecb1b4a

14
examples/mesh/tensionY_mono.load
View File

@ -1,14 +0,0 @@
# initial elastic step
$Loadcase 1 t 0.0005 N 1 f_out 1
Face 3 Y -0.025
Face 4 X 0.0
Face 4 Y 0.0
Face 4 Z 0.0
$EndLoadcase
# plastic step
$Loadcase 2 t 1.0 N 10 f_out 2
Face 3 Y -0.025
Face 4 X 0.0
Face 4 Y 0.0
Face 4 Z 0.0
$EndLoadcase

22
examples/mesh/tensionY_mono.yaml Normal file
View File

@ -0,0 +1,22 @@
---
loadstep:
- boundary_conditions:
mechanical:
- dot_u: ['x', -0.025, 'x']
tag: 3
- dot_u: [0.0, 0.0, 0.0]
tag: 4
discretization:
t: 0.0005
N: 1
f_out: 1
- boundary_conditions:
mechanical:
- dot_u: ['x', -0.025, 'x']
tag: 3
- dot_u: [0.0, 0.0, 0.0]
tag: 4
discretization:
t: 1.0
N: 10
f_out: 2

14
examples/mesh/tensionZ_3g.load
View File

@ -1,14 +0,0 @@
# initial elastic step
$Loadcase 1 t 0.0005 N 1 f_out 1
Face 1 Z 0.01
Face 2 X 0.0
Face 2 Y 0.0
Face 2 Z 0.0
$EndLoadcase
# plastic step
$Loadcase 2 t 1.0 N 10 f_out 2
Face 1 Z 0.01
Face 2 X 0.0
Face 2 Y 0.0
Face 2 Z 0.0
$EndLoadcase

22
examples/mesh/tensionZ_3g.yaml Normal file
View File

@ -0,0 +1,22 @@
---
loadstep:
- boundary_conditions:
mechanical:
- dot_u: ['x', 'x', 0.01]
tag: 1
- dot_u: [0.0, 0.0, 0.0]
tag: 2
discretization:
t: 0.0005
N: 1
f_out: 1
- boundary_conditions:
mechanical:
- dot_u: ['x', 'x', 0.01]
tag: 1
- dot_u: [0.0, 0.0, 0.0]
tag: 2
discretization:
t: 1.0
N: 10
f_out: 2

18
examples/mesh/tensionZ_5g.load
View File

@ -1,18 +0,0 @@
# initial elastic step
$Loadcase 1 t 0.0005 N 1 f_out 1
Face 1 X 0.0
Face 1 Y 0.0
Face 1 Z 0.0
Face 2 X 0.0
Face 2 Y 0.0
Face 2 Z 0.0025
$EndLoadcase
# plastic step
$Loadcase 2 t 1.0 N 10 f_out 2
Face 1 X 0.0
Face 1 Y 0.0
Face 1 Z 0.0
Face 2 X 0.0
Face 2 Y 0.0
Face 2 Z 0.0025
$EndLoadcase

22
examples/mesh/tensionZ_5g.yaml Normal file
View File

@ -0,0 +1,22 @@
---
loadstep:
- boundary_conditions:
mechanical:
- dot_u: [0.0, 0.0, 0.0]
tag: 1
- dot_u: [0.0, 0.0, 0.0025]
tag: 2
discretization:
t: 0.0005
N: 1
f_out: 1
- boundary_conditions:
mechanical:
- dot_u: [0.0, 0.0, 0.0]
tag: 1
- dot_u: [0.0, 0.0, 0.0025]
tag: 2
discretization:
t: 1.0
N: 10
f_out: 2

249
src/mesh/DAMASK_mesh.f90
View File

@ -23,38 +23,32 @@ program DAMASK_mesh
implicit none(type,external) implicit none(type,external)
type :: tLoadCase type :: tLoadCase
real(pREAL) :: time = 0.0_pREAL !< length of increment real(pREAL) :: t = 0.0_pREAL !< length of increment
integer :: incs = 0, & !< number of increments integer :: N = 0, & !< number of increments
outputfrequency = 1 !< frequency of result writes f_out = 1 !< frequency of result writes
logical :: followFormerTrajectory = .true. !< follow trajectory of former loadcase logical :: estimate_rate = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: faceID integer, allocatable, dimension(:) :: tag
type(tFieldBC), allocatable, dimension(:) :: fieldBC type(tMechBC) :: mechBC
end type tLoadCase end type tLoadCase
!--------------------------------------------------------------------------------------------------
! variables related to information from load case and geom file
integer, allocatable, dimension(:) :: chunkPos ! this is longer than needed for geometry parsing
integer :: &
N_def = 0 !< # of rate of deformation specifiers found in load case file
character(len=:), allocatable :: &
line
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! loop variables, convergence etc. ! loop variables, convergence etc.
integer, parameter :: & integer, parameter :: &
subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pREAL) :: & real(pREAL) :: &
time = 0.0_pREAL, & !< elapsed time t = 0.0_pREAL, & !< elapsed time
time0 = 0.0_pREAL, & !< begin of interval t_0 = 0.0_pREAL, & !< begin of interval
timeinc = 0.0_pREAL, & !< current time interval Delta_t = 0.0_pREAL, & !< current time interval
timeIncOld = 0.0_pREAL, & !< previous time interval Delta_t_prev = 0.0_pREAL, & !< previous time interval
remainingLoadCaseTime = 0.0_pREAL !< remaining time of current load case t_remaining = 0.0_pREAL !< remaining time of current load case
logical :: & logical :: &
guess, & !< guess along former trajectory guess, & !< guess along former trajectory
stagIterate stagIterate
integer :: & integer :: &
l, & l, &
i, & i, &
m, &
errorID, & errorID, &
cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$ cutBackLevel = 0, & !< cut back level \f$ t = \frac{t_{inc}}{2^l} \f$
stepFraction = 0, & !< fraction of current time interval stepFraction = 0, & !< fraction of current time interval
@ -67,8 +61,16 @@ program DAMASK_mesh
component component
type(tDict), pointer :: & type(tDict), pointer :: &
num_solver, & num_solver, &
num_mesh num_mesh, &
character(len=pSTRLEN), dimension(:), allocatable :: fileContent load, &
load_step, &
step_bc, &
mech_BC, &
step_discretization
type(tList), pointer :: &
load_steps, &
mech_u, &
step_mech
character(len=pSTRLEN) :: & character(len=pSTRLEN) :: &
incInfo, & incInfo, &
loadcase_string loadcase_string
@ -83,6 +85,9 @@ program DAMASK_mesh
integer(kind(COMPONENT_UNDEFINED_ID)) :: ID integer(kind(COMPONENT_UNDEFINED_ID)) :: ID
external :: & external :: &
quit quit
character(len=:), allocatable :: &
fileContent, fname
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! init DAMASK (all modules) ! init DAMASK (all modules)
@ -104,135 +109,101 @@ program DAMASK_mesh
CHKERRA(err_PETSc) CHKERRA(err_PETSc)
allocate(solres(1)) allocate(solres(1))
!-------------------------------------------------------------------------------------------------- if (worldrank == 0) then
! reading basic information from load case file and allocate data structure containing load cases fileContent = IO_read(CLI_loadFile)
fileContent = IO_readlines(trim(CLI_loadFile)) fname = CLI_loadFile
do l = 1, size(fileContent) if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
line = fileContent(l) call result_openJobFile(parallel=.false.)
if (IO_isBlank(line)) cycle ! skip empty lines call result_addSetupFile(fileContent,fname,'load case definition (mesh solver)')
call result_closeJobFile()
end if
chunkPos = IO_strPos(line) call parallelization_bcast_str(fileContent)
do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase load => YAML_parse_str_asDict(fileContent)
select case (IO_strValue(line,chunkPos,i)) load_steps => load%get_list('loadstep')
case('$Loadcase')
N_def = N_def + 1
end select
end do ! count all identifiers to allocate memory and do sanity check
end do
if (N_def < 1) call IO_error(error_ID = 837) allocate(loadCases(load_steps%length))
allocate(loadCases(N_def))
do i = 1, size(loadCases) do l = 1, load_steps%length
allocate(loadCases(i)%fieldBC(1)) load_step => load_steps%get_dict(l)
loadCases(i)%fieldBC(1)%ID = FIELD_MECH_ID step_bc => load_step%get_dict('boundary_conditions')
end do step_mech => step_bc%get_list('mechanical')
loadCases(l)%mechBC%nComponents = dimPlex !< X, Y (, Z) displacements
do i = 1, size(loadCases) allocate(loadCases(l)%mechBC%componentBC(dimPlex))
loadCases(i)%fieldBC(1)%nComponents = dimPlex !< X, Y (, Z) displacements do component = 1, dimPlex
allocate(loadCases(i)%fieldBC(1)%componentBC(loadCases(i)%fieldBC(1)%nComponents))
do component = 1, loadCases(i)%fieldBC(1)%nComponents
select case (component) select case (component)
case (1) case (1)
loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_X_ID loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_X_ID
case (2) case (2)
loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Y_ID loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_Y_ID
case (3) case (3)
loadCases(i)%fieldBC(1)%componentBC(component)%ID = COMPONENT_MECH_Z_ID loadCases(l)%mechBC%componentBC(component)%ID = COMPONENT_MECH_Z_ID
end select end select
end do end do
do component = 1, loadCases(i)%fieldBC(1)%nComponents do component = 1, dimPlex
allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pREAL) allocate(loadCases(l)%mechBC%componentBC(component)%Value(mesh_Nboundaries), source = 0.0_pREAL)
allocate(loadCases(i)%fieldBC(1)%componentBC(component)%Mask (mesh_Nboundaries), source = .false.) allocate(loadCases(l)%mechBC%componentBC(component)%Mask (mesh_Nboundaries), source = .false.)
end do end do
end do
!-------------------------------------------------------------------------------------------------- do m = 1, step_mech%length
! reading the load case and assign values to the allocated data structure mech_BC => step_mech%get_dict(m)
do l = 1, size(fileContent) currentFaceSet = -1
line = fileContent(l) do faceSet = 1, mesh_Nboundaries
if (IO_isBlank(line)) cycle ! skip empty lines if (mesh_boundaries(faceSet) == mech_BC%get_asInt('tag')) currentFaceSet = faceSet
end do
chunkPos = IO_strPos(line) if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
do i = 1, chunkPos(1) do component = 1, dimPlex
select case (IO_strValue(line,chunkPos,i)) mech_u => mech_BC%get_list('dot_u')
!-------------------------------------------------------------------------------------------------- if (mech_u%get_asStr(component) /= 'x') then
! loadcase information loadCases(l)%mechBC%componentBC(component)%Mask(currentFaceSet) = .true.
case('$Loadcase') loadCases(l)%mechBC%componentBC(component)%Value(currentFaceSet) = mech_u%get_asReal(component)
currentLoadCase = IO_intValue(line,chunkPos,i+1) end if
case('Face') end do
currentFace = IO_intValue(line,chunkPos,i+1)
currentFaceSet = -1
do faceSet = 1, mesh_Nboundaries
if (mesh_boundaries(faceSet) == currentFace) currentFaceSet = faceSet
end do
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
case('t')
loadCases(currentLoadCase)%time = IO_realValue(line,chunkPos,i+1)
case('N')
loadCases(currentLoadCase)%incs = IO_intValue(line,chunkPos,i+1)
case('f_out')
loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,chunkPos,i+1)
case('estimate_rate')
loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory
!--------------------------------------------------------------------------------------------------
! boundary condition information
case('X','Y','Z')
select case(IO_strValue(line,chunkPos,i))
case('X')
ID = COMPONENT_MECH_X_ID
case('Y')
ID = COMPONENT_MECH_Y_ID
case('Z')
ID = COMPONENT_MECH_Z_ID
end select
do component = 1, loadcases(currentLoadCase)%fieldBC(1)%nComponents
if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%ID == ID) then
loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask (currentFaceSet) = &
.true.
loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(currentFaceSet) = &
IO_realValue(line,chunkPos,i+1)
end if
end do
end select
end do end do
step_discretization => load_step%get_dict('discretization')
loadCases(l)%t = step_discretization%get_asReal('t')
loadCases(l)%N = step_discretization%get_asInt ('N')
if (load_step%get_asStr('f_out',defaultVal='n/a') == 'none') then
loadCases(l)%f_out = huge(0)
else
loadCases(l)%f_out = load_step%get_asInt('f_out', defaultVal=1)
end if
loadCases(l)%estimate_rate = (load_step%get_asBool('estimate_rate',defaultVal=.true.) .and. l>1)
end do end do
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! consistency checks and output of load case ! consistency checks and output of load case
loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase
errorID = 0 errorID = 0
checkLoadcases: do currentLoadCase = 1, size(loadCases) checkLoadcases: do l = 1, load_steps%length
write (loadcase_string, '(i0)' ) currentLoadCase write (loadcase_string, '(i0)' ) l
print'(/,1x,a,1x,i0)', 'load case:', currentLoadCase print'(/,1x,a,1x,i0)', 'load case:', l
if (.not. loadCases(currentLoadCase)%followFormerTrajectory) & if (.not. loadCases(l)%estimate_rate) &
print'(2x,a)', 'drop guessing along trajectory' print'(2x,a)', 'drop guessing along trajectory'
print'(2x,a)', 'Field '//trim(FIELD_MECH_label) print'(2x,a)', 'Field '//trim(FIELD_MECH_label)
do faceSet = 1, mesh_Nboundaries do faceSet = 1, mesh_Nboundaries
do component = 1, loadCases(currentLoadCase)%fieldBC(1)%nComponents do component = 1, loadCases(l)%mechBC%nComponents
if (loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Mask(faceSet)) & if (loadCases(l)%mechBC%componentBC(component)%Mask(faceSet)) &
print'(a,i2,a,i2,a,f12.7)', & print'(a,i2,a,i2,a,f12.7)', &
' Face ', mesh_boundaries(faceSet), & ' Face ', mesh_boundaries(faceSet), &
' Component ', component, & ' Component ', component, &
' Value ', loadCases(currentLoadCase)%fieldBC(1)%componentBC(component)%Value(faceSet) ' Value ', loadCases(l)%mechBC%componentBC(component)%Value(faceSet)
end do end do
end do end do
print'(2x,a,f12.6)', 'time: ', loadCases(currentLoadCase)%time print'(2x,a,f12.6)', 'time: ', loadCases(l)%t
if (loadCases(currentLoadCase)%incs < 1) errorID = 835 ! non-positive incs count if (loadCases(l)%N < 1) errorID = 835 ! non-positive incs count
print'(2x,a,i5)', 'increments: ', loadCases(currentLoadCase)%incs print'(2x,a,i5)', 'increments: ', loadCases(l)%N
if (loadCases(currentLoadCase)%outputfrequency < 1) errorID = 836 ! non-positive result frequency if (loadCases(l)%f_out < 1) errorID = 836 ! non-positive result frequency
print'(2x,a,i5)', 'output frequency: ', & print'(2x,a,i5)', 'output frequency: ', &
loadCases(currentLoadCase)%outputfrequency loadCases(l)%f_out
if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message if (errorID > 0) call IO_error(error_ID = errorID, ext_msg = loadcase_string) ! exit with error message
end do checkLoadcases end do checkLoadcases
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! doing initialization depending on active solvers ! doing initialization depending on active solvers
call FEM_Utilities_init(num_mesh) call FEM_Utilities_init(num_mesh)
call FEM_mechanical_init(loadCases(1)%fieldBC(1),num_mesh) call FEM_mechanical_init(loadCases(1)%mechBC,num_mesh)
call config_numerics_deallocate() call config_numerics_deallocate()
if (worldrank == 0) then if (worldrank == 0) then
@ -244,46 +215,46 @@ program DAMASK_mesh
flush(IO_STDOUT) flush(IO_STDOUT)
call materialpoint_result(0,0.0_pREAL) call materialpoint_result(0,0.0_pREAL)
loadCaseLooping: do currentLoadCase = 1, size(loadCases) loadCaseLooping: do l = 1, load_steps%length
time0 = time ! load case start time t_0 = t ! load case start time
guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1, loadCases(currentLoadCase)%incs incLooping: do inc = 1, loadCases(l)%N
totalIncsCounter = totalIncsCounter + 1 totalIncsCounter = totalIncsCounter + 1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end Delta_t_prev = Delta_t ! last timeinc that brought former inc to an end
timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pREAL) Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL)
timeinc = timeinc * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step Delta_t = Delta_t * real(subStepFactor,pREAL)**real(-cutBackLevel,pREAL) ! depending on cut back level, decrease time step
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time t_remaining = loadCases(l)%t + t_0 - t
time = time + timeinc ! forward target time t = t + Delta_t ! forward target time
stepFraction = stepFraction + 1 ! count step stepFraction = stepFraction + 1 ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
print'(/,1x,a)', '###########################################################################' print'(/,1x,a)', '###########################################################################'
print'(1x,a,es12.5,6(a,i0))',& print'(1x,a,es12.5,6(a,i0))',&
'Time', time, & 'Time', t, &
's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& 's: Increment ', inc, '/', loadCases(l)%N,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,& '-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', currentLoadCase,'/',size(loadCases) ' of load case ', l,'/',load_steps%length
write(incInfo,'(4(a,i0))') & write(incInfo,'(4(a,i0))') &
'Increment ',totalIncsCounter,'/',sum(loadCases%incs),& 'Increment ',totalIncsCounter,'/',sum(loadCases%N),&
'-',stepFraction, '/', subStepFactor**cutBackLevel '-',stepFraction, '/', subStepFactor**cutBackLevel
flush(IO_STDOUT) flush(IO_STDOUT)
call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) call FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,loadCases(l)%mechBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! solve fields ! solve fields
stagIter = 0 stagIter = 0
stagIterate = .true. stagIterate = .true.
do while (stagIterate) do while (stagIterate)
solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(currentLoadCase)%fieldBC(1)) solres(1) = FEM_mechanical_solution(incInfo,Delta_t,Delta_t_prev,loadCases(l)%mechBC)
if (.not. solres(1)%converged) exit if (.not. solres(1)%converged) exit
stagIter = stagIter + 1 stagIter = stagIter + 1
@ -294,13 +265,13 @@ program DAMASK_mesh
! check solution ! check solution
cutBack = .False. cutBack = .False.
if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found if (.not. all(solres(:)%converged .and. solres(:)%stagConverged)) then ! no solution found
if (cutBackLevel < maxCutBack) then ! do cut back if (cutBackLevel < maxCutBack) then ! do cut back
cutBack = .True. cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1 cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time t = t - Delta_t ! rewind time
timeinc = timeinc/2.0_pREAL Delta_t = Delta_t/2.0_pREAL
print'(/,1x,a)', 'cutting back' print'(/,1x,a)', 'cutting back'
else ! default behavior, exit if spectral solver does not converge else ! default behavior, exit if spectral solver does not converge
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
@ -308,10 +279,10 @@ program DAMASK_mesh
end if end if
else else
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc Delta_t_prev = Delta_t
end if end if
if (.not. cutBack .and. worldrank == 0) then if (.not. cutBack .and. worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, & write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
solres%converged, solres%iterationsNeeded ! write statistics about accepted solution solres%converged, solres%iterationsNeeded ! write statistics about accepted solution
flush(statUnit) flush(statUnit)
end if end if
@ -325,10 +296,10 @@ program DAMASK_mesh
print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged' print'(/,1x,a,1x,i0,1x,a)', 'increment', totalIncsCounter, 'NOT converged'
end if; flush(IO_STDOUT) end if; flush(IO_STDOUT)
if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency if (mod(inc,loadCases(l)%f_out) == 0) then ! at output frequency
print'(/,1x,a)', '... saving results ........................................................' print'(/,1x,a)', '... saving results ........................................................'
call FEM_mechanical_updateCoords() call FEM_mechanical_updateCoords()
call materialpoint_result(totalIncsCounter,time) call materialpoint_result(totalIncsCounter,t)
end if end if

24
src/mesh/FEM_utilities.f90
View File

@ -38,10 +38,6 @@ module FEM_utilities
character(len=*), parameter, public :: & character(len=*), parameter, public :: &
FIELD_MECH_label = 'mechanical' FIELD_MECH_label = 'mechanical'
enum, bind(c); enumerator :: &
FIELD_UNDEFINED_ID, &
FIELD_MECH_ID
end enum
enum, bind(c); enumerator :: & enum, bind(c); enumerator :: &
COMPONENT_UNDEFINED_ID, & COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, & COMPONENT_MECH_X_ID, &
@ -64,11 +60,10 @@ module FEM_utilities
logical, allocatable, dimension(:) :: Mask logical, allocatable, dimension(:) :: Mask
end type tComponentBC end type tComponentBC
type, public :: tFieldBC type, public :: tMechBC
integer(kind(FIELD_UNDEFINED_ID)) :: ID
integer :: nComponents = 0 integer :: nComponents = 0
type(tComponentBC), allocatable, dimension(:) :: componentBC type(tComponentBC), allocatable, dimension(:) :: componentBC
end type tFieldBC end type tMechBC
external :: & ! ToDo: write interfaces external :: & ! ToDo: write interfaces
PetscSectionGetFieldComponents, & PetscSectionGetFieldComponents, &
@ -79,7 +74,6 @@ module FEM_utilities
FEM_utilities_init, & FEM_utilities_init, &
utilities_constitutiveResponse, & utilities_constitutiveResponse, &
utilities_projectBCValues, & utilities_projectBCValues, &
FIELD_MECH_ID, &
COMPONENT_UNDEFINED_ID, & COMPONENT_UNDEFINED_ID, &
COMPONENT_MECH_X_ID, & COMPONENT_MECH_X_ID, &
COMPONENT_MECH_Y_ID, & COMPONENT_MECH_Y_ID, &
@ -142,9 +136,9 @@ end subroutine FEM_utilities_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates constitutive response !> @brief calculates constitutive response
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData) subroutine utilities_constitutiveResponse(Delta_t,P_av,forwardData)
real(pREAL), intent(in) :: timeinc !< loading time real(pREAL), intent(in) :: Delta_t !< loading time
logical, intent(in) :: forwardData !< age results logical, intent(in) :: forwardData !< age results
real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress
@ -152,9 +146,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
print'(/,1x,a)', '... evaluating constitutive response ......................................' print'(/,1x,a)', '... evaluating constitutive response ......................................'
call homogenization_mechanical_response(timeinc,1,mesh_maxNips*mesh_NcpElems) ! calculate P field call homogenization_mechanical_response(Delta_t,1,mesh_maxNips*mesh_NcpElems) ! calculate P field
if (.not. terminallyIll) & if (.not. terminallyIll) &
call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems]) call homogenization_mechanical_response2(Delta_t,[1,mesh_maxNips],[1,mesh_NcpElems])
cutBack = .false. cutBack = .false.
P_av = sum(homogenization_P,dim=3) * wgt P_av = sum(homogenization_P,dim=3) * wgt
@ -168,7 +162,7 @@ end subroutine utilities_constitutiveResponse
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Project BC values to local vector !> @brief Project BC values to local vector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,timeinc) subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCValue,BCDotValue,Delta_t)
Vec :: localVec Vec :: localVec
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
@ -176,7 +170,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
IS :: bcPointsIS IS :: bcPointsIS
PetscInt, pointer :: bcPoints(:) PetscInt, pointer :: bcPoints(:)
real(pREAL), pointer :: localArray(:) real(pREAL), pointer :: localArray(:)
real(pREAL) :: BCValue,BCDotValue,timeinc real(pREAL) :: BCValue,BCDotValue,Delta_t
PetscErrorCode :: err_PETSc PetscErrorCode :: err_PETSc
@ -193,7 +187,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,err_PETSc) call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do dof = offset+comp+1, offset+numDof, numComp do dof = offset+comp+1, offset+numDof, numComp
localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc localArray(dof) = localArray(dof) + BCValue + BCDotValue*Delta_t
end do end do
end do end do
call VecRestoreArrayF90(localVec,localArray,err_PETSc) call VecRestoreArrayF90(localVec,localArray,err_PETSc)

54
src/mesh/mesh_mech_FEM.f90
View File

@ -36,8 +36,8 @@ module mesh_mechanical_FEM
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! derived types ! derived types
type tSolutionParams type tSolutionParams
type(tFieldBC) :: fieldBC type(tMechBC) :: mechBC
real(pREAL) :: timeinc real(pREAL) :: Delta_t
end type tSolutionParams end type tSolutionParams
type(tSolutionParams) :: params type(tSolutionParams) :: params
@ -97,9 +97,9 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields and fills them with data !> @brief allocates all neccessary fields and fills them with data
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_init(fieldBC,num_mesh) subroutine FEM_mechanical_init(mechBC,num_mesh)
type(tFieldBC), intent(in) :: fieldBC type(tMechBC), intent(in) :: mechBC
type(tDict), pointer, intent(in) :: num_mesh type(tDict), pointer, intent(in) :: num_mesh
DM :: mechanical_mesh DM :: mechanical_mesh
@ -209,14 +209,14 @@ subroutine FEM_mechanical_init(fieldBC,num_mesh)
end do end do
numBC = 0 numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1 if (mechBC%componentBC(field)%Mask(faceSet)) numBC = numBC + 1
end do; end do end do; end do
allocate(pbcField(numBC), source=0_pPETSCINT) allocate(pbcField(numBC), source=0_pPETSCINT)
allocate(pbcComps(numBC)) allocate(pbcComps(numBC))
allocate(pbcPoints(numBC)) allocate(pbcPoints(numBC))
numBC = 0 numBC = 0
do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries do field = 1, dimPlex; do faceSet = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(faceSet)) then if (mechBC%componentBC(field)%Mask(faceSet)) then
numBC = numBC + 1 numBC = numBC + 1
call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc) call ISCreateGeneral(PETSC_COMM_WORLD,1_pPETSCINT,[field-1],PETSC_COPY_VALUES,pbcComps(numBC),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -320,15 +320,15 @@ end subroutine FEM_mechanical_init
!> @brief solution for the FEM load step !> @brief solution for the FEM load step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
type(tSolutionState) function FEM_mechanical_solution( & type(tSolutionState) function FEM_mechanical_solution( &
incInfoIn,timeinc,timeinc_old,fieldBC) incInfoIn,Delta_t,Delta_t_prev,mechBC)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! input data for solution ! input data for solution
real(pREAL), intent(in) :: & real(pREAL), intent(in) :: &
timeinc, & !< increment in time for current solution Delta_t, & !< increment in time for current solution
timeinc_old !< increment in time of last increment Delta_t_prev !< increment in time of last increment
type(tFieldBC), intent(in) :: & type(tMechBC), intent(in) :: &
fieldBC mechBC
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
incInfoIn incInfoIn
@ -339,8 +339,8 @@ type(tSolutionState) function FEM_mechanical_solution( &
FEM_mechanical_solution%converged = .false. FEM_mechanical_solution%converged = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%Delta_t = Delta_t
params%fieldBC = fieldBC params%mechBC = mechBC
call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,err_PETSc) ! solve mechanical_snes based on solution guess (result in solution) call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,err_PETSc) ! solve mechanical_snes based on solution guess (result in solution)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
@ -407,13 +407,13 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries do field = 1_pPETSCINT, dimPlex; do face = 1_pPETSCINT, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then if (params%mechBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc) call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end if end if
@ -459,7 +459,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(params%timeinc,P_av,ForwardData) call utilities_constitutiveResponse(params%Delta_t,P_av,ForwardData)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,err_MPI) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1_MPI_INTEGER_KIND,MPI_LOGICAL,MPI_LOR,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'
ForwardData = .false. ForwardData = .false.
@ -557,13 +557,13 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc) call VecWAXPY(x_local,1.0_pREAL,xx_local,solution_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (params%fieldBC%componentBC(field)%Mask(face)) then if (params%mechBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, & call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,params%fieldBC%componentBC(field)%Value(face),params%timeinc) 0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc) call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end if end if
@ -665,13 +665,13 @@ end subroutine FEM_mechanical_formJacobian
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief forwarding routine !> @brief forwarding routine
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC) subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
type(tFieldBC), intent(in) :: & type(tMechBC), intent(in) :: &
fieldBC mechBC
real(pREAL), intent(in) :: & real(pREAL), intent(in) :: &
timeinc_old, & Delta_t_prev, &
timeinc Delta_t
logical, intent(in) :: & logical, intent(in) :: &
guess guess
@ -702,13 +702,13 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
call VecAXPY(solution_local,1.0_pREAL,x_local,err_PETSc) call VecAXPY(solution_local,1.0_pREAL,x_local,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
do field = 1, dimPlex; do face = 1, mesh_Nboundaries do field = 1, dimPlex; do face = 1, mesh_Nboundaries
if (fieldBC%componentBC(field)%Mask(face)) then if (mechBC%componentBC(field)%Mask(face)) then
call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc) call DMGetStratumSize(dm_local,'Face Sets',mesh_boundaries(face),bcSize,err_PETSc)
if (bcSize > 0) then if (bcSize > 0) then
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc) call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, & call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,fieldBC%componentBC(field)%Value(face),timeinc_old) 0.0_pREAL,mechBC%componentBC(field)%Value(face),Delta_t_prev)
call ISDestroy(bcPoints,err_PETSc) call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end if end if
@ -721,12 +721,12 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,fieldBC)
! update rate and forward last inc ! update rate and forward last inc
call VecCopy(solution,solution_rate,err_PETSc) call VecCopy(solution,solution_rate,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecScale(solution_rate,timeinc_old**(-1),err_PETSc) call VecScale(solution_rate,Delta_t_prev**(-1),err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end if end if
call VecCopy(solution_rate,solution,err_PETSc) call VecCopy(solution_rate,solution,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecScale(solution,timeinc,err_PETSc) call VecScale(solution,Delta_t,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
end subroutine FEM_mechanical_forward end subroutine FEM_mechanical_forward