better names

This commit is contained in:
Sharan Roongta 2023-12-16 18:23:21 +01:00
parent 8828783bd4
commit e4ad5a94a5
4 changed files with 47 additions and 54 deletions

@ -1 +1 @@
Subproject commit cc7b00fe7b1c195e1b68a2a9c0e77c4b533af847
Subproject commit 117b65d852659158c0f4ca3bf8ce8db51a7a1961

View File

@ -23,32 +23,25 @@ program DAMASK_mesh
implicit none(type,external)
type :: tLoadCase
real(pREAL) :: t = 0.0_pREAL !< length of increment
integer :: N = 0, & !< number of increments
f_out = 1 !< frequency of result writes
logical :: estimate_rate = .true. !< follow trajectory of former loadcase
real(pREAL) :: t = 0.0_pREAL !< length of increment
integer :: N = 0, & !< number of increments
f_out = 1 !< frequency of result writes
logical :: estimate_rate = .true. !< follow trajectory of former loadcase
integer, allocatable, dimension(:) :: tag
type(tMechBC) :: mechBC
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.
integer, parameter :: &
subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pREAL) :: &
time = 0.0_pREAL, & !< elapsed time
time0 = 0.0_pREAL, & !< begin of interval
timeinc = 0.0_pREAL, & !< current time interval
timeIncOld = 0.0_pREAL, & !< previous time interval
remainingLoadCaseTime = 0.0_pREAL !< remaining time of current load case
t = 0.0_pREAL, & !< elapsed time
t_0 = 0.0_pREAL, & !< begin of interval
Delta_t = 0.0_pREAL, & !< current time interval
Delta_t_prev = 0.0_pREAL, & !< previous time interval
t_remaining = 0.0_pREAL !< remaining time of current load case
logical :: &
guess, & !< guess along former trajectory
stagIterate
@ -160,7 +153,7 @@ program DAMASK_mesh
end do
if (currentFaceSet < 0) call IO_error(error_ID = 837, ext_msg = 'invalid BC')
do component = 1, dimPlex
mech_u => mech_BC%get_list('u')
mech_u => mech_BC%get_list('dot_u')
if (mech_u%get_asStr(component) /= 'x') then
loadCases(l)%mechBC%componentBC(component)%Mask(currentFaceSet) = .true.
loadCases(l)%mechBC%componentBC(component)%Value(currentFaceSet) = mech_u%get_asReal(component)
@ -223,7 +216,7 @@ program DAMASK_mesh
call materialpoint_result(0,0.0_pREAL)
loadCaseLooping: do l = 1, load_steps%length
time0 = time ! load case start time
t_0 = t ! load case start time
guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1, loadCases(l)%N
@ -231,21 +224,21 @@ program DAMASK_mesh
!--------------------------------------------------------------------------------------------------
! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end
timeinc = 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_prev = Delta_t ! last timeinc that brought former inc to an end
Delta_t = loadCases(l)%t/real(loadCases(l)%N,pREAL)
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
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(l)%t + time0 - time
time = time + timeinc ! forward target time
t_remaining = loadCases(l)%t + t_0 - t
t = t + Delta_t ! forward target time
stepFraction = stepFraction + 1 ! count step
!--------------------------------------------------------------------------------------------------
! report begin of new step
print'(/,1x,a)', '###########################################################################'
print'(1x,a,es12.5,6(a,i0))',&
'Time', time, &
'Time', t, &
's: Increment ', inc, '/', loadCases(l)%N,&
'-', stepFraction, '/', subStepFactor**cutBackLevel,&
' of load case ', l,'/',load_steps%length
@ -254,14 +247,14 @@ program DAMASK_mesh
'-',stepFraction, '/', subStepFactor**cutBackLevel
flush(IO_STDOUT)
call FEM_mechanical_forward(guess,timeinc,timeIncOld,loadCases(l)%mechBC)
call FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,loadCases(l)%mechBC)
!--------------------------------------------------------------------------------------------------
! solve fields
stagIter = 0
stagIterate = .true.
do while (stagIterate)
solres(1) = FEM_mechanical_solution(incInfo,timeinc,timeIncOld,loadCases(l)%mechBC)
solres(1) = FEM_mechanical_solution(incInfo,Delta_t,Delta_t_prev,loadCases(l)%mechBC)
if (.not. solres(1)%converged) exit
stagIter = stagIter + 1
@ -272,13 +265,13 @@ program DAMASK_mesh
! check solution
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
cutBack = .True.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time
timeinc = timeinc/2.0_pREAL
t = t - Delta_t ! rewind time
Delta_t = Delta_t/2.0_pREAL
print'(/,1x,a)', 'cutting back'
else ! default behavior, exit if spectral solver does not converge
if (worldrank == 0) close(statUnit)
@ -286,10 +279,10 @@ program DAMASK_mesh
end if
else
guess = .true. ! start guessing after first converged (sub)inc
timeIncOld = timeinc
Delta_t_prev = Delta_t
end if
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
flush(statUnit)
end if
@ -306,7 +299,7 @@ program DAMASK_mesh
if (mod(inc,loadCases(l)%f_out) == 0) then ! at output frequency
print'(/,1x,a)', '... saving results ........................................................'
call FEM_mechanical_updateCoords()
call materialpoint_result(totalIncsCounter,time)
call materialpoint_result(totalIncsCounter,t)
end if

View File

@ -136,9 +136,9 @@ end subroutine FEM_utilities_init
!--------------------------------------------------------------------------------------------------
!> @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
real(pREAL),intent(out), dimension(3,3) :: P_av !< average PK stress
@ -146,9 +146,9 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
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) &
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.
P_av = sum(homogenization_P,dim=3) * wgt
@ -162,7 +162,7 @@ end subroutine utilities_constitutiveResponse
!--------------------------------------------------------------------------------------------------
!> @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
PetscInt :: field, comp, nBcPoints, point, dof, numDof, numComp, offset
@ -170,7 +170,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
IS :: bcPointsIS
PetscInt, pointer :: bcPoints(:)
real(pREAL), pointer :: localArray(:)
real(pREAL) :: BCValue,BCDotValue,timeinc
real(pREAL) :: BCValue,BCDotValue,Delta_t
PetscErrorCode :: err_PETSc
@ -187,7 +187,7 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
call PetscSectionGetFieldOffset(section,bcPoints(point),field,offset,err_PETSc)
CHKERRQ(err_PETSc)
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
call VecRestoreArrayF90(localVec,localArray,err_PETSc)

View File

@ -37,7 +37,7 @@ module mesh_mechanical_FEM
! derived types
type tSolutionParams
type(tMechBC) :: mechBC
real(pREAL) :: timeinc
real(pREAL) :: Delta_t
end type tSolutionParams
type(tSolutionParams) :: params
@ -320,13 +320,13 @@ end subroutine FEM_mechanical_init
!> @brief solution for the FEM load step
!--------------------------------------------------------------------------------------------------
type(tSolutionState) function FEM_mechanical_solution( &
incInfoIn,timeinc,timeinc_old,mechBC)
incInfoIn,Delta_t,Delta_t_prev,mechBC)
!--------------------------------------------------------------------------------------------------
! input data for solution
real(pREAL), intent(in) :: &
timeinc, & !< increment in time for current solution
timeinc_old !< increment in time of last increment
Delta_t, & !< increment in time for current solution
Delta_t_prev !< increment in time of last increment
type(tMechBC), intent(in) :: &
mechBC
character(len=*), intent(in) :: &
@ -339,7 +339,7 @@ type(tSolutionState) function FEM_mechanical_solution( &
FEM_mechanical_solution%converged = .false.
!--------------------------------------------------------------------------------------------------
! set module wide availabe data
params%timeinc = timeinc
params%Delta_t = Delta_t
params%mechBC = mechBC
call SNESSolve(mechanical_snes,PETSC_NULL_VEC,solution,err_PETSc) ! solve mechanical_snes based on solution guess (result in solution)
@ -413,7 +413,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%timeinc)
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
end if
@ -459,7 +459,7 @@ subroutine FEM_mechanical_formResidual(dm_local,xx_local,f_local,dummy,err_PETSc
!--------------------------------------------------------------------------------------------------
! 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)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
ForwardData = .false.
@ -563,7 +563,7 @@ subroutine FEM_mechanical_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,err_P
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_projectBCValues(x_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%timeinc)
0.0_pREAL,params%mechBC%componentBC(field)%Value(face),params%Delta_t)
call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
end if
@ -665,13 +665,13 @@ end subroutine FEM_mechanical_formJacobian
!--------------------------------------------------------------------------------------------------
!> @brief forwarding routine
!--------------------------------------------------------------------------------------------------
subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,mechBC)
subroutine FEM_mechanical_forward(guess,Delta_t,Delta_t_prev,mechBC)
type(tMechBC), intent(in) :: &
mechBC
real(pREAL), intent(in) :: &
timeinc_old, &
timeinc
Delta_t_prev, &
Delta_t
logical, intent(in) :: &
guess
@ -708,7 +708,7 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,mechBC)
call DMGetStratumIS(dm_local,'Face Sets',mesh_boundaries(face),bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
call utilities_projectBCValues(solution_local,section,0_pPETSCINT,field-1,bcPoints, &
0.0_pREAL,mechBC%componentBC(field)%Value(face),timeinc_old)
0.0_pREAL,mechBC%componentBC(field)%Value(face),Delta_t_prev)
call ISDestroy(bcPoints,err_PETSc)
CHKERRQ(err_PETSc)
end if
@ -721,12 +721,12 @@ subroutine FEM_mechanical_forward(guess,timeinc,timeinc_old,mechBC)
! update rate and forward last inc
call VecCopy(solution,solution_rate,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)
end if
call VecCopy(solution_rate,solution,err_PETSc)
CHKERRQ(err_PETSc)
call VecScale(solution,timeinc,err_PETSc)
call VecScale(solution,Delta_t,err_PETSc)
CHKERRQ(err_PETSc)
end subroutine FEM_mechanical_forward