From a2e942033607a9d764aa5d64154001b364d7b5ff Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Sep 2020 14:58:44 +0200 Subject: [PATCH 01/15] boundary conditions to not change during iteration --- src/grid/DAMASK_grid.f90 | 90 ++++++++++---------- src/grid/grid_mech_FEM.f90 | 21 +++-- src/grid/grid_mech_spectral_basic.f90 | 21 ++--- src/grid/grid_mech_spectral_polarisation.f90 | 21 ++--- 4 files changed, 73 insertions(+), 80 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 0e30025d8..08f2e3f88 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -23,7 +23,7 @@ program DAMASK_grid use grid_damage_spectral use grid_thermal_spectral use results - + implicit none !-------------------------------------------------------------------------------------------------- @@ -87,7 +87,7 @@ program DAMASK_grid mech_updateCoords procedure(grid_mech_spectral_basic_restartWrite), pointer :: & mech_restartWrite - + external :: & quit class (tNode), pointer :: & @@ -96,10 +96,10 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) - - call CPFEM_initAll + + call CPFEM_initAll write(6,'(/,a)') ' <<<+- DAMASK_spectral init -+>>>'; flush(6) - + write(6,'(/,a)') ' Shanthraj et al., Handbook of Mechanics of Materials, 2019' write(6,'(a)') ' https://doi.org/10.1007/978-981-10-6855-3_80' @@ -122,16 +122,16 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type - + debug_grid => debug_root%get('grid',defaultVal=emptyList) - select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) + select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) case ('Basic') mech_init => grid_mech_spectral_basic_init mech_forward => grid_mech_spectral_basic_forward mech_solution => grid_mech_spectral_basic_solution mech_updateCoords => grid_mech_spectral_basic_updateCoords mech_restartWrite => grid_mech_spectral_basic_restartWrite - + case ('Polarisation') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -140,7 +140,7 @@ program DAMASK_grid mech_solution => grid_mech_spectral_polarisation_solution mech_updateCoords => grid_mech_spectral_polarisation_updateCoords mech_restartWrite => grid_mech_spectral_polarisation_restartWrite - + case ('FEM') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -149,24 +149,24 @@ program DAMASK_grid mech_solution => grid_mech_FEM_solution mech_updateCoords => grid_mech_FEM_updateCoords mech_restartWrite => grid_mech_FEM_restartWrite - + case default call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) - + end select - + !-------------------------------------------------------------------------------------------------- -! reading information from load case file and to sanity checks +! reading information from load case file and to sanity checks fileContent = IO_readlines(trim(loadCaseFile)) if(size(fileContent) == 0) call IO_error(307,ext_msg='No load case specified') - + allocate (loadCases(0)) ! array of load cases do currentLoadCase = 1, size(fileContent) line = fileContent(currentLoadCase) if (IO_isBlank(line)) cycle chunkPos = IO_stringPos(line) - + do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase select case (IO_lc(IO_stringValue(line,chunkPos,i))) case('l','fdot','dotf','f') @@ -179,7 +179,7 @@ program DAMASK_grid enddo if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase - + newLoadCase%stress%myType='stress' field = 1 newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default @@ -191,7 +191,7 @@ program DAMASK_grid field = field + 1 newLoadCase%ID(field) = FIELD_DAMAGE_ID endif damageActive - + call newLoadCase%rot%fromEulers(real([0.0,0.0,0.0],pReal)) readIn: do i = 1, chunkPos(1) select case (IO_lc(IO_stringValue(line,chunkPos,i))) @@ -257,9 +257,9 @@ program DAMASK_grid call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector)) end select enddo readIn - + newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1) ! by default, guess from previous load case - + reportAndCheck: if (worldrank == 0) then write (loadcase_string, '(i0)' ) currentLoadCase write(6,'(/,1x,a,i0)') 'load case: ', currentLoadCase @@ -324,13 +324,13 @@ program DAMASK_grid select case (loadCases(1)%ID(field)) case(FIELD_MECH_ID) call mech_init - + case(FIELD_THERMAL_ID) call grid_thermal_spectral_init - + case(FIELD_DAMAGE_ID) call grid_damage_spectral_init - + end select enddo @@ -348,16 +348,16 @@ program DAMASK_grid '.sta',form='FORMATTED', position='APPEND', status='OLD') endif writeHeader endif - + writeUndeformed: if (interface_restartInc < 1) then write(6,'(1/,a)') ' ... writing initial configuration to file ........................' call CPFEM_results(0,0.0_pReal) endif writeUndeformed - + loadCaseLooping: do currentLoadCase = 1, size(loadCases) time0 = time ! load case start time guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc - + incLooping: do inc = 1, loadCases(currentLoadCase)%incs totalIncsCounter = totalIncsCounter + 1 @@ -382,13 +382,13 @@ program DAMASK_grid endif endif timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step - + skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc? time = time + timeinc ! just advance time, skip already performed calculation guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference else skipping stepFraction = 0 ! fraction scaled by stepFactor**cutLevel - + subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time time = time + timeinc ! forward target time @@ -417,7 +417,7 @@ program DAMASK_grid deformation_BC = loadCases(currentLoadCase)%deformation, & stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rot) - + case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) end select @@ -432,21 +432,21 @@ program DAMASK_grid do field = 1, nActiveFields select case(loadCases(currentLoadCase)%ID(field)) case(FIELD_MECH_ID) - solres(field) = mech_solution (& - incInfo,timeinc,timeIncOld, & + solres(field) = mech_solution(& + incInfo, & stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rot) - + case(FIELD_THERMAL_ID) solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) - + case(FIELD_DAMAGE_ID) solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) - + end select - + if (.not. solres(field)%converged) exit ! no solution found - + enddo stagIter = stagIter + 1 stagIterate = stagIter < stagItMax & @@ -480,17 +480,17 @@ program DAMASK_grid if (worldrank == 0) close(statUnit) call quit(0) ! quit endif - + enddo subStepLooping - + cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc - + if (all(solres(:)%converged)) then write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' converged' else write(6,'(/,a,i0,a)') ' increment ', totalIncsCounter, ' NOT converged' endif; flush(6) - + if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency write(6,'(1/,a)') ' ... writing results to file ......................................' flush(6) @@ -501,17 +501,17 @@ program DAMASK_grid call CPFEM_restartWrite endif endif skipping - + enddo incLooping - + enddo loadCaseLooping - - + + !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation write(6,'(/,a)') ' ###########################################################################' if (worldrank == 0) close(statUnit) - + call quit(0) ! no complains ;) - + end program DAMASK_grid diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index e8e1345ef..f89e4fec9 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -268,15 +268,12 @@ end subroutine grid_mech_FEM_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the FEM scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_FEM_solution(incInfoIn,stress_BC,rotation_BC) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC type(rotation), intent(in) :: & @@ -293,13 +290,6 @@ function grid_mech_FEM_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -341,6 +331,15 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, PetscScalar, pointer, dimension(:,:,:,:) :: & u_current,u_lastInc +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + + call DMDAVecGetArrayF90(mech_grid,solution_current,u_current,ierr); CHKERRQ(ierr) call DMDAVecGetArrayF90(mech_grid,solution_lastInc,u_lastInc,ierr); CHKERRQ(ierr) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 66694e516..2b5e55e45 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -232,15 +232,12 @@ end subroutine grid_mech_spectral_basic_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the basic scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_basic_solution(incInfoIn,stress_BC,rotation_BC) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC type(rotation), intent(in) :: & @@ -259,14 +256,6 @@ function grid_mech_spectral_basic_solution(incInfoIn,timeinc,timeinc_old,stress_ S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old - !-------------------------------------------------------------------------------------------------- ! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -306,6 +295,14 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo PetscErrorCode :: ierr PetscScalar, dimension(:,:,:,:), pointer :: F +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) if (cutBack) then diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index f6d1b9ebb..0cd40545d 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -259,15 +259,12 @@ end subroutine grid_mech_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_polarisation_solution(incInfoIn,stress_BC,rotation_BC) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - real(pReal), intent(in) :: & - timeinc, & !< time increment of current solution - timeinc_old !< time increment of last successful increment type(tBoundaryCondition), intent(in) :: & stress_BC type(rotation), intent(in) :: & @@ -290,14 +287,6 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,timeinc,timeinc_old, S_scale = math_invSym3333(C_minMaxAvg) endif -!-------------------------------------------------------------------------------------------------- -! set module wide available data - params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values - params%rotation_BC = rotation_BC - params%timeinc = timeinc - params%timeincOld = timeinc_old - !-------------------------------------------------------------------------------------------------- ! solve BVP call SNESsolve(snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) @@ -339,6 +328,14 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc integer :: i, j, k real(pReal), dimension(3,3) :: F_lambda33 +!-------------------------------------------------------------------------------------------------- +! set module wide available data + params%stress_mask = stress_BC%maskFloat + params%stress_BC = stress_BC%values + params%rotation_BC = rotation_BC + params%timeinc = timeinc + params%timeincOld = timeinc_old + call DMDAVecGetArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) F => FandF_tau(0: 8,:,:,:) F_tau => FandF_tau(9:17,:,:,:) From 053c3f39ea1e0bd014baa20c93cf70934e7b1fd0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 19 Sep 2020 23:40:17 +0200 Subject: [PATCH 02/15] solution completely relies on state defined by 'forward' --- src/grid/DAMASK_grid.f90 | 8 +------- src/grid/grid_mech_FEM.f90 | 8 ++------ src/grid/grid_mech_spectral_basic.f90 | 8 ++------ src/grid/grid_mech_spectral_polarisation.f90 | 8 ++------ 4 files changed, 7 insertions(+), 25 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 878ed868b..7bc737ecd 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -433,17 +433,11 @@ program DAMASK_grid do field = 1, nActiveFields select case(loadCases(currentLoadCase)%ID(field)) case(FIELD_MECH_ID) - solres(field) = mech_solution (& - incInfo, & - stress_BC = loadCases(currentLoadCase)%stress, & - rotation_BC = loadCases(currentLoadCase)%rot) - + solres(field) = mech_solution(incInfo) case(FIELD_THERMAL_ID) solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) - case(FIELD_DAMAGE_ID) solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) - end select if (.not. solres(field)%converged) exit ! no solution found diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 749c1b43e..18f66e9c0 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -269,16 +269,12 @@ end subroutine grid_mech_FEM_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the FEM scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_FEM_solution(incInfoIn,stress_BC,rotation_BC) result(solution) +function grid_mech_FEM_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -290,7 +286,7 @@ function grid_mech_FEM_solution(incInfoIn,stress_BC,rotation_BC) result(solution !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) !-------------------------------------------------------------------------------------------------- ! solve BVP diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 33ee1ff97..3b814ee44 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -233,16 +233,12 @@ end subroutine grid_mech_spectral_basic_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the basic scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_basic_solution(incInfoIn,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_basic_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -254,7 +250,7 @@ function grid_mech_spectral_basic_solution(incInfoIn,stress_BC,rotation_BC) resu !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 5470cf0b3..3c02d62c0 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -260,16 +260,12 @@ end subroutine grid_mech_spectral_polarisation_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Polarisation scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_mech_spectral_polarisation_solution(incInfoIn,stress_BC,rotation_BC) result(solution) +function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! input data for solution character(len=*), intent(in) :: & incInfoIn - type(tBoundaryCondition), intent(in) :: & - stress_BC - type(rotation), intent(in) :: & - rotation_BC type(tSolutionState) :: & solution !-------------------------------------------------------------------------------------------------- @@ -281,7 +277,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn,stress_BC,rotation_B !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(rotation_BC,stress_BC%maskLogical,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) if (num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg From 8dfb972ac11e4392e2f07090402e339e6903c994 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 11:49:20 +0200 Subject: [PATCH 03/15] private is already default (module wide) --- src/grid/grid_mech_FEM.f90 | 36 ++++++++++---------- src/grid/grid_mech_spectral_basic.f90 | 6 ++-- src/grid/grid_mech_spectral_polarisation.f90 | 2 +- src/grid/spectral_utilities.f90 | 16 ++++----- 4 files changed, 30 insertions(+), 30 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 4e181cf09..57afecb20 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -27,9 +27,9 @@ module grid_mech_FEM !-------------------------------------------------------------------------------------------------- ! derived types - type(tSolutionParams), private :: params + type(tSolutionParams) :: params - type, private :: tNumerics + type :: tNumerics integer :: & itmin, & !< minimum number of iterations itmax !< maximum number of iterations @@ -43,45 +43,45 @@ module grid_mech_FEM eps_stress_rtol !< relative tolerance for fullfillment of stress BC end type tNumerics - type(tNumerics), private :: num - logical, private:: & + type(tNumerics) :: num + logical :: & debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: mech_grid - SNES, private :: mech_snes - Vec, private :: solution_current, solution_lastInc, solution_rate + DM :: mech_grid + SNES :: mech_snes + Vec :: solution_current, solution_lastInc, solution_rate !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc - real(pReal), private :: detJ - real(pReal), private, dimension(3) :: delta - real(pReal), private, dimension(3,8) :: BMat - real(pReal), private, dimension(8,8) :: HGMat - PetscInt, private :: xstart,ystart,zstart,xend,yend,zend + real(pReal), dimension(:,:,:,:,:), allocatable :: F, P_current, F_lastInc + real(pReal) :: detJ + real(pReal), dimension(3) :: delta + real(pReal), dimension(3,8) :: BMat + real(pReal), dimension(8,8) :: HGMat + PetscInt :: xstart,ystart,zstart,xend,yend,zend !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & + real(pReal), dimension(3,3) :: & F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastIter = math_I3, & F_aim_lastInc = math_I3, & !< previous average deformation gradient P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - character(len=:), allocatable, private :: incInfo !< time and increment information + character(len=:), allocatable :: incInfo !< time and increment information - real(pReal), private, dimension(3,3,3,3) :: & + real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness S = 0.0_pReal !< current compliance (filled up with zeros) - real(pReal), private :: & + real(pReal) :: & err_BC !< deviation from stress BC - integer, private :: & + integer :: & totalIter = 0 !< total iteration in current increment public :: & diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 133b8f3ed..8ec6172b4 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -28,7 +28,7 @@ module grid_mech_spectral_basic ! derived types type(tSolutionParams) :: params - type, private :: tNumerics + type :: tNumerics logical :: update_gamma !< update gamma operator with current stiffness integer :: & itmin, & !< minimum number of iterations @@ -42,7 +42,7 @@ module grid_mech_spectral_basic type(tNumerics) :: num ! numerics parameters. Better name? - logical, private :: debugRotation + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -65,7 +65,7 @@ module grid_mech_spectral_basic P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress character(len=:), allocatable :: incInfo !< time and increment information - real(pReal), private, dimension(3,3,3,3) :: & + real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness C_minMaxAvg = 0.0_pReal, & !< current (min+max)/2 stiffness diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 475f278a4..c61f71b27 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -48,7 +48,7 @@ module grid_mech_spectral_polarisation type(tNumerics) :: num ! numerics parameters. Better name? - logical, private :: debugRotation + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 5e53fe2a5..52df30c39 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -47,15 +47,15 @@ module spectral_utilities complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:,:), pointer :: vectorField_fourier !< vector field fourier representation for fftw real(C_DOUBLE), public, dimension(:,:,:), pointer :: scalarField_real !< scalar field real representation for fftw complex(C_DOUBLE_COMPLEX),public, dimension(:,:,:), pointer :: scalarField_fourier !< scalar field fourier representation for fftw - complex(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method - complex(pReal), private, dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives - complex(pReal), private, dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives - real(pReal), private, dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness + complex(pReal), dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat !< gamma operator (field) for spectral method + complex(pReal), dimension(:,:,:,:), allocatable :: xi1st !< wave vector field for first derivatives + complex(pReal), dimension(:,:,:,:), allocatable :: xi2nd !< wave vector field for second derivatives + real(pReal), dimension(3,3,3,3) :: C_ref !< mechanic reference stiffness !-------------------------------------------------------------------------------------------------- ! plans for FFTW - type(C_PTR), private :: & + type(C_PTR) :: & planTensorForth, & !< FFTW MPI plan P(x) to P(k) planTensorBack, & !< FFTW MPI plan F(k) to F(x) planVectorForth, & !< FFTW MPI plan v(x) to v(k) @@ -65,7 +65,7 @@ module spectral_utilities !-------------------------------------------------------------------------------------------------- ! variables controlling debugging - logical, private :: & + 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 @@ -108,14 +108,14 @@ module spectral_utilities real(pReal) :: timeincOld end type tSolutionParams - type, private :: tNumerics + type :: tNumerics integer :: & divergence_correction !< scale divergence/curl calculation: [0: no correction, 1: size scaled to 1, 2: size scaled to Npoints] logical :: & memory_efficient !< calculate gamma operator on the fly end type tNumerics - type(tNumerics), private :: num ! numerics parameters. Better name? + type(tNumerics) :: num ! numerics parameters. Better name? enum, bind(c); enumerator :: & DERIVATIVE_CONTINUOUS_ID, & From d584207e0abc29156ec2fe35b529a5bd9adb7942 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 12:11:48 +0200 Subject: [PATCH 04/15] same layout for easy diff --- src/grid/grid_mech_FEM.f90 | 23 ++++++------- src/grid/grid_mech_spectral_basic.f90 | 34 ++++++++++---------- src/grid/grid_mech_spectral_polarisation.f90 | 25 +++++++------- 3 files changed, 39 insertions(+), 43 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 57afecb20..bd6f862b2 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -34,18 +34,15 @@ module grid_mech_FEM itmin, & !< minimum number of iterations itmax !< maximum number of iterations real(pReal) :: & - err_div, & - divTol, & - BCTol, & eps_div_atol, & !< absolute tolerance for equilibrium eps_div_rtol, & !< relative tolerance for equilibrium eps_stress_atol, & !< absolute tolerance for fullfillment of stress BC eps_stress_rtol !< relative tolerance for fullfillment of stress BC end type tNumerics - type(tNumerics) :: num - logical :: & - debugRotation + type(tNumerics) :: num ! numerics parameters. Better name? + + logical :: debugRotation !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -72,7 +69,6 @@ module grid_mech_FEM P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress character(len=:), allocatable :: incInfo !< time and increment information - real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness C_volAvgLastInc = 0.0_pReal, & !< previous volume average stiffness @@ -99,7 +95,6 @@ contains subroutine grid_mech_FEM_init real(pReal) :: HGCoeff = 0.0e-2_pReal - PetscInt, dimension(0:worldsize-1) :: localK real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal real(pReal), dimension(4,8) :: & @@ -111,24 +106,25 @@ subroutine grid_mech_FEM_init -1.0_pReal, 1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal,-1.0_pReal,-1.0_pReal,-1.0_pReal, & 1.0_pReal, 1.0_pReal, 1.0_pReal, 1.0_pReal], [4,8]) + real(pReal), dimension(3,3,3,3) :: devNull PetscErrorCode :: ierr + PetscScalar, pointer, dimension(:,:,:,:) :: & + u_current,u_lastInc + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle character(len=pStringLen) :: & fileName class(tNode), pointer :: & num_grid, & debug_grid - real(pReal), dimension(3,3,3,3) :: devNull - PetscScalar, pointer, dimension(:,:,:,:) :: & - u_current,u_lastInc print'(/,a)', ' <<<+- grid_mech_FEM init -+>>>'; flush(6) -!----------------------------------------------------------------------------------------------- +!------------------------------------------------------------------------------------------------- ! debugging options debug_grid => config_debug%get('grid', defaultVal=emptyList) debugRotation = debug_grid%contains('rotation') - + !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) @@ -242,6 +238,7 @@ subroutine grid_mech_FEM_init F_lastInc = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) ! initialize to identity F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) endif restartRead + materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent call utilities_updateCoords(F) call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 8ec6172b4..b199c723f 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -96,18 +96,17 @@ subroutine grid_mech_spectral_basic_init real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - class (tNode), pointer :: & - num_grid, & - debug_grid - PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & F ! pointer to solution data - PetscInt, dimension(worldsize) :: localK + PetscInt, dimension(0:worldsize-1) :: localK integer(HID_T) :: fileHandle, groupHandle integer :: fileUnit character(len=pStringLen) :: & fileName + class (tNode), pointer :: & + num_grid, & + debug_grid print'(/,a)', ' <<<+- grid_mech_spectral_basic init -+>>>'; flush(6) @@ -121,7 +120,7 @@ subroutine grid_mech_spectral_basic_init ! debugging options debug_grid => config_debug%get('grid', defaultVal=emptyList) debugRotation = debug_grid%contains('rotation') - + !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) @@ -140,7 +139,7 @@ subroutine grid_mech_spectral_basic_init if (num%eps_stress_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_stress_rtol') if (num%itmax <= 1) call IO_error(301,ext_msg='itmax') if (num%itmin > num%itmax .or. num%itmin < 1) call IO_error(301,ext_msg='itmin') - + !-------------------------------------------------------------------------------------------------- ! set default and user defined options for PETSc call PETScOptionsInsertString(PETSC_NULL_OPTIONS,'-mech_snes_type ngmres',ierr) @@ -157,8 +156,8 @@ subroutine grid_mech_spectral_basic_init ! initialize solver specific parts of PETSc call SNESCreate(PETSC_COMM_WORLD,snes,ierr); CHKERRQ(ierr) call SNESSetOptionsPrefix(snes,'mech_',ierr);CHKERRQ(ierr) - localK = 0 - localK(worldrank+1) = grid3 + localK = 0 + localK(worldrank) = grid3 call MPI_Allreduce(MPI_IN_PLACE,localK,worldsize,MPI_INTEGER,MPI_SUM,PETSC_COMM_WORLD,ierr) call DMDACreate3d(PETSC_COMM_WORLD, & DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, & ! cut off stencil at boundary @@ -202,8 +201,8 @@ subroutine grid_mech_spectral_basic_init endif restartRead materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateCoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_updateCoords(reshape(F,shape(F_lastInc))) + call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) ! deassociate pointer @@ -288,7 +287,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: F + PetscScalar, pointer, dimension(:,:,:,:) :: F !-------------------------------------------------------------------------------------------------- ! set module wide available data @@ -334,7 +333,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc - F = reshape(Utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -368,7 +367,7 @@ subroutine grid_mech_spectral_basic_restartWrite call DMDAVecGetArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) - print'(a)', ' writing solver data required for restart to file'; flush(6) + print*, 'writing solver data required for restart to file'; flush(6) write(fileName,'(a,a,i0,a)') trim(getSolverJobName()),'_',worldrank,'.hdf5' fileHandle = HDF5_openFile(fileName,'w') @@ -462,6 +461,7 @@ subroutine formResidual(in, F, & call SNESGetIterationNumber(snes,PETScIter,ierr); CHKERRQ(ierr) if (nfuncs == 0 .and. PETScIter == 0) totalIter = -1 ! new increment + !-------------------------------------------------------------------------------------------------- ! begin of new iteration newIteration: if (totalIter <= PETScIter) then @@ -484,8 +484,8 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim + deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) ! S = 0.0 for no bc + F_aim = F_aim -deltaF_aim err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- @@ -493,7 +493,7 @@ subroutine formResidual(in, F, & tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residuum ! store fPK field for subsequent FFT forward transform call utilities_FFTtensorForward ! FFT forward of global "tensorField_real" - err_div = Utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use + err_div = utilities_divergenceRMS() ! divRMS of tensorField_fourier for later use call utilities_fourierGammaConvolution(params%rotation_BC%rotate(deltaF_aim,active=.true.)) ! convolution of Gamma and tensorField_fourier call utilities_FFTtensorBackward ! FFT backward of global tensorField_fourier diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index c61f71b27..393f8f3f7 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -15,7 +15,6 @@ module grid_mech_spectral_polarisation use DAMASK_interface use HDF5_utilities use math - use rotations use spectral_utilities use FEsolving use config @@ -108,10 +107,6 @@ subroutine grid_mech_spectral_polarisation_init real(pReal), dimension(3,3,grid(1),grid(2),grid3) :: P real(pReal), dimension(3,3) :: & temp33_Real = 0.0_pReal - class (tNode), pointer :: & - num_grid, & - debug_grid - PetscErrorCode :: ierr PetscScalar, pointer, dimension(:,:,:,:) :: & FandF_tau, & ! overall pointer to solution data @@ -122,13 +117,16 @@ subroutine grid_mech_spectral_polarisation_init integer :: fileUnit character(len=pStringLen) :: & fileName + class (tNode), pointer :: & + num_grid, & + debug_grid print'(/,a)', ' <<<+- grid_mech_spectral_polarisation init -+>>>'; flush(6) print*, 'Shanthraj et al., International Journal of Plasticity 66:31–45, 2015' print*, 'https://doi.org/10.1016/j.ijplas.2014.02.006' -!------------------------------------------------------------------------------------------------ +!------------------------------------------------------------------------------------------------- ! debugging options debug_grid => config_debug%get('grid',defaultVal=emptyList) debugRotation = debug_grid%contains('rotation') @@ -136,6 +134,7 @@ subroutine grid_mech_spectral_polarisation_init !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) @@ -228,8 +227,8 @@ subroutine grid_mech_spectral_polarisation_init endif restartRead materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent - call Utilities_updateCoords(reshape(F,shape(F_lastInc))) - call Utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 + call utilities_updateCoords(reshape(F,shape(F_lastInc))) + call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 reshape(F,shape(F_lastInc)), & ! target F 0.0_pReal) ! time increment call DMDAVecRestoreArrayF90(da,solution_vec,FandF_tau,ierr); CHKERRQ(ierr) ! deassociate pointer @@ -277,7 +276,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) - if (num%update_gamma) then + if(num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg S_scale = math_invSym3333(C_minMaxAvg) @@ -320,7 +319,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc type(rotation), intent(in) :: & rotation_BC PetscErrorCode :: ierr - PetscScalar, dimension(:,:,:,:), pointer :: FandF_tau, F, F_tau + PetscScalar, pointer, dimension(:,:,:,:) :: FandF_tau, F, F_tau integer :: i, j, k real(pReal), dimension(3,3) :: F_lambda33 @@ -588,7 +587,7 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - params%stress_BC))) ! S = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, P_av - params%stress_BC) ! S = 0.0 for no bc err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & -params%rotation_BC%rotate(F_av)) + & params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc @@ -596,7 +595,7 @@ subroutine formResidual(in, FandF_tau, & tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise call utilities_FFTtensorForward - err_div = Utilities_divergenceRMS() !< root mean squared error in divergence of stress + err_div = utilities_divergenceRMS() !< root mean squared error in divergence of stress !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -615,7 +614,7 @@ subroutine formResidual(in, FandF_tau, & tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = F call utilities_FFTtensorForward - err_curl = Utilities_curlRMS() + err_curl = utilities_curlRMS() end subroutine formResidual From 6367cb8fcba017330d762273e4ee38fc514fe187 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 13:29:51 +0200 Subject: [PATCH 05/15] consistent handling of boundary conditions A stress boundary condition 'P' indicates the stress at the end of the load case (same as for 'F') 'Pdot' for given increase of (technical) strain is not implemented. Does not change anything for the most common case of zero-stress boundary conditions, but simplifies the specification of stress ramps --- src/grid/DAMASK_grid.f90 | 6 +++--- src/grid/grid_mech_FEM.f90 | 19 +++++++++--------- src/grid/grid_mech_spectral_basic.f90 | 21 ++++++++++++-------- src/grid/grid_mech_spectral_polarisation.f90 | 15 +++++++++----- 4 files changed, 36 insertions(+), 25 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 59a5b452c..54d080744 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -107,8 +107,8 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! initialize field solver information nActiveFields = 1 - if (any(thermal_type == THERMAL_conduction_ID )) nActiveFields = nActiveFields + 1 - if (any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 + if(any(thermal_type == THERMAL_conduction_ID)) nActiveFields = nActiveFields + 1 + if(any(damage_type == DAMAGE_nonlocal_ID )) nActiveFields = nActiveFields + 1 allocate(solres(nActiveFields)) allocate(newLoadCase%ID(nActiveFields)) @@ -181,7 +181,7 @@ program DAMASK_grid if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase - newLoadCase%stress%myType='stress' + newLoadCase%stress%myType='p' field = 1 newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default thermalActive: if (any(thermal_type == THERMAL_conduction_ID)) then diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index bd6f862b2..9efd88505 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -66,8 +66,8 @@ module grid_mech_FEM F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastIter = math_I3, & F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -327,7 +327,6 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -374,6 +373,12 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + endif + call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) call DMDAVecRestoreArrayF90(mech_grid,solution_current,u_current,ierr);CHKERRQ(ierr) @@ -489,8 +494,6 @@ subroutine formResidual(da_local,x_local, & PetscScalar, pointer,dimension(:,:,:,:) :: x_scal, f_scal PetscScalar, dimension(8,3) :: x_elem, f_elem PetscInt :: i, ii, j, jj, k, kk, ctr, ele - real(pReal), dimension(3,3) :: & - deltaF_aim PetscInt :: & PETScIter, & nfuncs @@ -539,10 +542,8 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim_lastIter = F_aim - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) - F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index b199c723f..cd795f005 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -62,8 +62,8 @@ module grid_mech_spectral_basic F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -129,7 +129,7 @@ subroutine grid_mech_spectral_basic_init num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.01_pReal) + num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.001_pReal) num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) @@ -292,7 +292,6 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -333,6 +332,12 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + endif + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) call DMDAVecRestoreArrayF90(da,solution_vec,F,ierr); CHKERRQ(ierr) @@ -484,9 +489,9 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - deltaF_aim = math_mul3333xx33(S, P_av - params%stress_BC) ! S = 0.0 for no bc - F_aim = F_aim -deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - params%stress_BC))) ! mask = 0.0 when no stress bc + deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc + F_aim = F_aim - deltaF_aim + err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme @@ -499,7 +504,7 @@ subroutine formResidual(in, F, & !-------------------------------------------------------------------------------------------------- ! constructing residual - residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too + residuum = tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) ! Gamma*P gives correction towards div(P) = 0, so needs to be zero, too end subroutine formResidual diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 393f8f3f7..4dc9d6d7c 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -70,8 +70,8 @@ module grid_mech_spectral_polarisation F_aim = math_I3, & !< current prescribed deformation gradient F_aim_lastInc = math_I3, & !< previous average deformation gradient F_av = 0.0_pReal, & !< average incompatible def grad field - P_av = 0.0_pReal !< average 1st Piola--Kirchhoff stress - + P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress + P_aim = 0.0_pReal character(len=:), allocatable :: incInfo !< time and increment information real(pReal), dimension(3,3,3,3) :: & C_volAvg = 0.0_pReal, & !< current volume average stiffness @@ -326,7 +326,6 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! set module wide available data params%stress_mask = stress_BC%maskFloat - params%stress_BC = stress_BC%values params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -373,6 +372,12 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc + if (stress_BC%myType=='p') then + P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + elseif (stress_BC%myType=='pdot') then !UNTESTED + P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + endif + F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average rotation_BC%rotate(F_aim,active=.true.)),& [9,grid(1),grid(2),grid3]) @@ -587,10 +592,10 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, P_av - params%stress_BC) ! S = 0.0 for no bc + F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & -params%rotation_BC%rotate(F_av)) + & - params%stress_mask * (P_av-params%stress_BC))) ! mask = 0.0 for no bc + params%stress_mask * (P_av-P_aim))) ! mask = 0.0 for no bc ! calculate divergence tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise From 0a7d4f61acfce1de1c7deb2b7962802641b0163b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 16:24:08 +0200 Subject: [PATCH 06/15] Need only logical mask 'merge' substitutes multiplication with float mask --- src/grid/DAMASK_grid.f90 | 105 +++++++++---------- src/grid/grid_mech_FEM.f90 | 26 +++-- src/grid/grid_mech_spectral_basic.f90 | 26 +++-- src/grid/grid_mech_spectral_polarisation.f90 | 30 +++--- src/grid/spectral_utilities.f90 | 15 ++- 5 files changed, 96 insertions(+), 106 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index 54d080744..d33d136f7 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -24,7 +24,7 @@ program DAMASK_grid use grid_damage_spectral use grid_thermal_spectral use results - + implicit none !-------------------------------------------------------------------------------------------------- @@ -88,7 +88,7 @@ program DAMASK_grid mech_updateCoords procedure(grid_mech_spectral_basic_restartWrite), pointer :: & mech_restartWrite - + external :: & quit class (tNode), pointer :: & @@ -97,10 +97,10 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) - - call CPFEM_initAll + + call CPFEM_initAll print'(/,a)', ' <<<+- DAMASK_spectral init -+>>>'; flush(6) - + print*, 'Shanthraj et al., Handbook of Mechanics of Materials, 2019' print*, 'https://doi.org/10.1007/978-981-10-6855-3_80' @@ -123,16 +123,16 @@ program DAMASK_grid !-------------------------------------------------------------------------------------------------- ! assign mechanics solver depending on selected type - + debug_grid => config_debug%get('grid',defaultVal=emptyList) - select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) + select case (trim(num_grid%get_asString('solver', defaultVal = 'Basic'))) case ('Basic') mech_init => grid_mech_spectral_basic_init mech_forward => grid_mech_spectral_basic_forward mech_solution => grid_mech_spectral_basic_solution mech_updateCoords => grid_mech_spectral_basic_updateCoords mech_restartWrite => grid_mech_spectral_basic_restartWrite - + case ('Polarisation') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -141,7 +141,7 @@ program DAMASK_grid mech_solution => grid_mech_spectral_polarisation_solution mech_updateCoords => grid_mech_spectral_polarisation_updateCoords mech_restartWrite => grid_mech_spectral_polarisation_restartWrite - + case ('FEM') if(debug_grid%contains('basic')) & call IO_warning(42, ext_msg='debug Divergence') @@ -150,24 +150,24 @@ program DAMASK_grid mech_solution => grid_mech_FEM_solution mech_updateCoords => grid_mech_FEM_updateCoords mech_restartWrite => grid_mech_FEM_restartWrite - + case default call IO_error(error_ID = 891, ext_msg = trim(num_grid%get_asString('solver'))) - + end select - + !-------------------------------------------------------------------------------------------------- -! reading information from load case file and to sanity checks +! reading information from load case file and to sanity checks fileContent = IO_readlines(trim(interface_loadFile)) if(size(fileContent) == 0) call IO_error(307,ext_msg='No load case specified') - + allocate (loadCases(0)) ! array of load cases do currentLoadCase = 1, size(fileContent) line = fileContent(currentLoadCase) if (IO_isBlank(line)) cycle chunkPos = IO_stringPos(line) - + do i = 1, chunkPos(1) ! reading compulsory parameters for loadcase select case (IO_lc(IO_stringValue(line,chunkPos,i))) case('l','fdot','dotf','f') @@ -180,7 +180,7 @@ program DAMASK_grid enddo if ((N_def /= N_n) .or. (N_n /= N_t) .or. N_n < 1) & ! sanity check call IO_error(error_ID=837,el=currentLoadCase,ext_msg = trim(interface_loadFile)) ! error message for incomplete loadcase - + newLoadCase%stress%myType='p' field = 1 newLoadCase%ID(field) = FIELD_MECH_ID ! mechanical active by default @@ -192,7 +192,7 @@ program DAMASK_grid field = field + 1 newLoadCase%ID(field) = FIELD_DAMAGE_ID endif damageActive - + call newLoadCase%rot%fromEulers(real([0.0,0.0,0.0],pReal)) readIn: do i = 1, chunkPos(1) select case (IO_lc(IO_stringValue(line,chunkPos,i))) @@ -210,18 +210,16 @@ program DAMASK_grid temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not a * if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - newLoadCase%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) ! logical mask in 3x3 notation - newLoadCase%deformation%maskFloat = merge(ones,zeros,newLoadCase%deformation%maskLogical) ! float (1.0/0.0) mask in 3x3 notation - newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation + newLoadCase%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation + newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation case('p','stress', 's') temp_valueVector = 0.0_pReal do j = 1, 9 temp_maskVector(j) = IO_stringValue(line,chunkPos,i+j) /= '*' ! true if not an asterisk if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo - newLoadCase%stress%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) - newLoadCase%stress%maskFloat = merge(ones,zeros,newLoadCase%stress%maskLogical) - newLoadCase%stress%values = math_9to33(temp_valueVector) + newLoadCase%stress%mask = transpose(reshape(temp_maskVector,[ 3,3])) + newLoadCase%stress%values = math_9to33(temp_valueVector) case('t','time','delta') ! increment time newLoadCase%time = IO_floatValue(line,chunkPos,i+1) case('n','incs','increments') ! number of increments @@ -258,9 +256,9 @@ program DAMASK_grid call newLoadCase%rot%fromMatrix(math_9to33(temp_valueVector)) end select enddo readIn - + newLoadCase%followFormerTrajectory = merge(.true.,.false.,currentLoadCase > 1) ! by default, guess from previous load case - + reportAndCheck: if (worldrank == 0) then write (loadcase_string, '(i0)' ) currentLoadCase print'(/,a,i0)', ' load case: ', currentLoadCase @@ -268,8 +266,8 @@ program DAMASK_grid print*, ' drop guessing along trajectory' if (newLoadCase%deformation%myType == 'l') then do j = 1, 3 - if (any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & - any(newLoadCase%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined + if (any(newLoadCase%deformation%mask(j,1:3) .eqv. .true.) .and. & + any(newLoadCase%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined enddo print*, ' velocity gradient:' else if (newLoadCase%deformation%myType == 'f') then @@ -278,20 +276,19 @@ program DAMASK_grid print*, ' deformation gradient rate:' endif do i = 1, 3; do j = 1, 3 - if(newLoadCase%deformation%maskLogical(i,j)) then + if(newLoadCase%deformation%mask(i,j)) then write(6,'(2x,f12.7)',advance='no') newLoadCase%deformation%values(i,j) else write(6,'(2x,12a)',advance='no') ' * ' endif enddo; write(6,'(/)',advance='no') enddo - if (any(newLoadCase%stress%maskLogical .eqv. & - newLoadCase%deformation%maskLogical)) errorID = 831 ! exclusive or masking only - if (any(newLoadCase%stress%maskLogical .and. transpose(newLoadCase%stress%maskLogical) & - .and. (math_I3<1))) errorID = 838 ! no rotation is allowed by stress BC + if (any(newLoadCase%stress%mask .eqv. newLoadCase%deformation%mask)) errorID = 831 ! exclusive or masking only + if (any(newLoadCase%stress%mask .and. transpose(newLoadCase%stress%mask) .and. (math_I3<1))) & + errorID = 838 ! no rotation is allowed by stress BC print*, ' stress / GPa:' do i = 1, 3; do j = 1, 3 - if(newLoadCase%stress%maskLogical(i,j)) then + if(newLoadCase%stress%mask(i,j)) then write(6,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal else write(6,'(2x,12a)',advance='no') ' * ' @@ -325,13 +322,13 @@ program DAMASK_grid select case (loadCases(1)%ID(field)) case(FIELD_MECH_ID) call mech_init - + case(FIELD_THERMAL_ID) call grid_thermal_spectral_init - + case(FIELD_DAMAGE_ID) call grid_damage_spectral_init - + end select enddo @@ -348,16 +345,16 @@ program DAMASK_grid '.sta',form='FORMATTED', position='APPEND', status='OLD') endif writeHeader endif - + writeUndeformed: if (interface_restartInc < 1) then print'(/,a)', ' ... writing initial configuration to file ........................' call CPFEM_results(0,0.0_pReal) endif writeUndeformed - + loadCaseLooping: do currentLoadCase = 1, size(loadCases) time0 = time ! load case start time guess = loadCases(currentLoadCase)%followFormerTrajectory ! change of load case? homogeneous guess for the first inc - + incLooping: do inc = 1, loadCases(currentLoadCase)%incs totalIncsCounter = totalIncsCounter + 1 @@ -382,13 +379,13 @@ program DAMASK_grid endif endif timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step - + skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc? time = time + timeinc ! just advance time, skip already performed calculation guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference else skipping stepFraction = 0 ! fraction scaled by stepFactor**cutLevel - + subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) remainingLoadCaseTime = loadCases(currentLoadCase)%time+time0 - time time = time + timeinc ! forward target time @@ -417,7 +414,7 @@ program DAMASK_grid deformation_BC = loadCases(currentLoadCase)%deformation, & stress_BC = loadCases(currentLoadCase)%stress, & rotation_BC = loadCases(currentLoadCase)%rot) - + case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) end select @@ -438,9 +435,9 @@ program DAMASK_grid case(FIELD_DAMAGE_ID) solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) end select - + if (.not. solres(field)%converged) exit ! no solution found - + enddo stagIter = stagIter + 1 stagIterate = stagIter < stagItMax & @@ -474,17 +471,17 @@ program DAMASK_grid if (worldrank == 0) close(statUnit) call quit(0) ! quit endif - + enddo subStepLooping - + cutBackLevel = max(0, cutBackLevel - 1) ! try half number of subincs next inc - + if (all(solres(:)%converged)) then print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' converged' else print'(/,a,i0,a)', ' increment ', totalIncsCounter, ' NOT converged' endif; flush(6) - + if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0) then ! at output frequency print'(1/,a)', ' ... writing results to file ......................................' flush(6) @@ -495,17 +492,17 @@ program DAMASK_grid call CPFEM_restartWrite endif endif skipping - + enddo incLooping - + enddo loadCaseLooping - - + + !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation print'(/,a)', ' ###########################################################################' if (worldrank == 0) close(statUnit) - + call quit(0) ! no complains ;) - + end program DAMASK_grid diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 9efd88505..9d9c962e1 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -25,8 +25,6 @@ module grid_mech_FEM implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -282,7 +280,7 @@ function grid_mech_FEM_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) !-------------------------------------------------------------------------------------------------- ! solve BVP @@ -326,7 +324,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat + params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -340,20 +338,20 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, else C_volAvgLastInc = C_volAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif if (guess) then @@ -374,9 +372,9 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime, ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc if (stress_BC%myType=='p') then - P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) elseif (stress_BC%myType=='pdot') then !UNTESTED - P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif call VecAXPY(solution_current,timeinc,solution_rate,ierr); CHKERRQ(ierr) @@ -543,7 +541,7 @@ subroutine formResidual(da_local,x_local, & !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc - err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! constructing residual diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index cd795f005..fd5e63663 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -24,8 +24,6 @@ module grid_mech_spectral_basic implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -247,7 +245,7 @@ function grid_mech_spectral_basic_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) if(num%update_gamma) call utilities_updateGamma(C_minMaxAvg) !-------------------------------------------------------------------------------------------------- @@ -291,7 +289,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat + params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -305,20 +303,20 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -333,9 +331,9 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc if (stress_BC%myType=='p') then - P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) elseif (stress_BC%myType=='pdot') then !UNTESTED - P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average @@ -491,7 +489,7 @@ subroutine formResidual(in, F, & ! stress BC handling deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - deltaF_aim - err_BC = maxval(abs(params%stress_mask * (P_av - P_aim))) ! mask = 0.0 when no stress bc + err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 4dc9d6d7c..91632f3b9 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -24,8 +24,6 @@ module grid_mech_spectral_polarisation implicit none private -!-------------------------------------------------------------------------------------------------- -! derived types type(tSolutionParams) :: params type :: tNumerics @@ -275,7 +273,7 @@ function grid_mech_spectral_polarisation_solution(incInfoIn) result(solution) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) - S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask>.5_pReal,C_volAvg) + S = utilities_maskedCompliance(params%rotation_BC,params%stress_mask,C_volAvg) if(num%update_gamma) then call utilities_updateGamma(C_minMaxAvg) C_scale = C_minMaxAvg @@ -325,7 +323,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc !-------------------------------------------------------------------------------------------------- ! set module wide available data - params%stress_mask = stress_BC%maskFloat + params%stress_mask = stress_BC%mask params%rotation_BC = rotation_BC params%timeinc = timeinc params%timeincOld = timeinc_old @@ -341,20 +339,20 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc C_volAvgLastInc = C_volAvg C_minMaxAvgLastInc = C_minMaxAvg - F_aimDot = merge(stress_BC%maskFloat*(F_aim-F_aim_lastInc)/timeinc_old, 0.0_pReal, guess) + F_aimDot = merge(merge((F_aim-F_aim_lastInc)/timeinc_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) F_aim_lastInc = F_aim !----------------------------------------------------------------------------------------------- ! calculate rate for aim if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc) + F_aimDot = F_aimDot & + + merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * deformation_BC%values + F_aimDot = F_aimDot & + + merge(deformation_BC%values,.0_pReal,deformation_BC%mask) elseif (deformation_BC%myType=='f') then ! aim at end of load case is prescribed - F_aimDot = & - F_aimDot + deformation_BC%maskFloat * (deformation_BC%values - F_aim_lastInc)/loadCaseTime + F_aimDot = F_aimDot & + + merge((deformation_BC%values - F_aim_lastInc)/loadCaseTime,.0_pReal,deformation_BC%mask) endif Fdot = utilities_calculateRate(guess, & @@ -373,9 +371,9 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc ! update average and local deformation gradients F_aim = F_aim_lastInc + F_aimDot * timeinc if (stress_BC%myType=='p') then - P_aim = P_aim + stress_BC%maskFloat*(stress_BC%values - P_aim)/loadCaseTime*timeinc + P_aim = P_aim + merge((stress_BC%values - P_aim)/loadCaseTime*timeinc,.0_pReal,stress_BC%mask) elseif (stress_BC%myType=='pdot') then !UNTESTED - P_aim = P_aim + stress_BC%maskFloat*stress_BC%values*timeinc + P_aim = P_aim + merge(stress_BC%values*timeinc,.0_pReal,stress_BC%mask) endif F = reshape(utilities_forwardField(timeinc,F_lastInc,Fdot, & ! estimate of F at end of time+timeinc that matches rotated F_aim on average @@ -593,9 +591,9 @@ subroutine formResidual(in, FandF_tau, & !-------------------------------------------------------------------------------------------------- ! stress BC handling F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc - err_BC = maxval(abs((1.0_pReal-params%stress_mask) * math_mul3333xx33(C_scale,F_aim & - -params%rotation_BC%rotate(F_av)) + & - params%stress_mask * (P_av-P_aim))) ! mask = 0.0 for no bc + err_BC = maxval(abs(merge(P_av-P_aim, & + math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& + params%stress_mask))) ! calculate divergence tensorField_real = 0.0_pReal tensorField_real(1:3,1:3,1:grid(1),1:grid(2),1:grid3) = residual_F !< stress field in disguise diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 52df30c39..896337bf6 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -82,10 +82,9 @@ module spectral_utilities end type tSolutionState type, public :: tBoundaryCondition !< set of parameters defining a boundary condition - real(pReal), dimension(3,3) :: values = 0.0_pReal, & - maskFloat = 0.0_pReal - logical, dimension(3,3) :: maskLogical = .false. - character(len=pStringLen) :: myType = 'None' + real(pReal), dimension(3,3) :: values = 0.0_pReal + logical, dimension(3,3) :: mask = .false. + character(len=pStringLen) :: myType = 'None' end type tBoundaryCondition type, public :: tLoadCase @@ -101,11 +100,11 @@ module spectral_utilities integer(kind(FIELD_UNDEFINED_ID)), allocatable :: ID(:) end type tLoadCase - type, public :: tSolutionParams !< @todo use here the type definition for a full loadcase - real(pReal), dimension(3,3) :: stress_mask, stress_BC + type, public :: tSolutionParams + real(pReal), dimension(3,3) :: stress_BC + logical, dimension(3,3) :: stress_mask type(rotation) :: rotation_BC - real(pReal) :: timeinc - real(pReal) :: timeincOld + real(pReal) :: timeinc, timeincOld end type tSolutionParams type :: tNumerics From c4d0ac71a999c8d10220f3748d84f35505fd6936 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 17:15:39 +0200 Subject: [PATCH 07/15] silences GNU10 compilation warning only for 'DEFENSIVE'. Does not cause any harm --- src/rotations.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index 72046a965..04490fe3c 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -1452,7 +1452,7 @@ subroutine selfTest contains - function quaternion_equal(qu1,qu2) result(ok) + pure recursive function quaternion_equal(qu1,qu2) result(ok) real(pReal), intent(in), dimension(4) :: qu1,qu2 logical :: ok From 329cc1c953cd24cccaa0002152f8cfdc546415e4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 17:36:11 +0200 Subject: [PATCH 08/15] tighter tolerance for stress --- src/grid/grid_mech_FEM.f90 | 12 +++++------ src/grid/grid_mech_spectral_basic.f90 | 14 ++++++------- src/grid/grid_mech_spectral_polarisation.f90 | 22 ++++++++++---------- 3 files changed, 24 insertions(+), 24 deletions(-) diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 9d9c962e1..a6f713adc 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -126,12 +126,12 @@ subroutine grid_mech_FEM_init !------------------------------------------------------------------------------------------------- ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) - num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) - num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) - num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) - num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal) + num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index fd5e63663..4ed41cb7b 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -123,13 +123,13 @@ subroutine grid_mech_spectral_basic_init ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) - num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) - num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) - num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=0.001_pReal) - num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) - num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) + num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal) + num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 91632f3b9..b8bbcacd8 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -133,17 +133,17 @@ subroutine grid_mech_spectral_polarisation_init ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) - num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) - num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) - num%eps_curl_atol = num_grid%get_asFloat ('eps_curl_atol', defaultVal=1.0e-10_pReal) - num%eps_curl_rtol = num_grid%get_asFloat ('eps_curl_rtol', defaultVal=5.0e-4_pReal) - num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol', defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol', defaultVal=0.01_pReal) - num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) - num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) - num%alpha = num_grid%get_asFloat ('alpha', defaultVal=1.0_pReal) - num%beta = num_grid%get_asFloat ('beta', defaultVal=1.0_pReal) + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) + num%eps_div_atol = num_grid%get_asFloat('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_curl_atol = num_grid%get_asFloat('eps_curl_atol', defaultVal=1.0e-10_pReal) + num%eps_curl_rtol = num_grid%get_asFloat('eps_curl_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal) + num%itmin = num_grid%get_asInt ('itmin', defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax', defaultVal=250) + num%alpha = num_grid%get_asFloat('alpha', defaultVal=1.0_pReal) + num%beta = num_grid%get_asFloat('beta', defaultVal=1.0_pReal) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') From 9ed92781833fdd25b6b95c23fa60f9c37bb828a2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 20 Sep 2020 19:46:33 +0200 Subject: [PATCH 09/15] polishing --- src/material.f90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index f8e487b66..cabc57835 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -2,7 +2,7 @@ !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Parses material config file, either solverJobName.materialConfig or material.config +!> @brief Defines phase and homogenization !-------------------------------------------------------------------------------------------------- module material use prec @@ -83,7 +83,7 @@ module material material_homogenizationAt !< homogenization ID of each element integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack material_homogenizationMemberAt !< position of the element within its homogenization instance - integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) + integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) material_phaseAt !< phase ID of each element integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,elem) material_phaseMemberAt !< position of the element within its phase instance @@ -326,7 +326,7 @@ subroutine material_parseMicrostructure constituents, & !> list of constituents constituent, & !> constituent definition phases, & - homogenization + homogenizations integer, dimension(:), allocatable :: & counterPhase, & @@ -362,8 +362,8 @@ subroutine material_parseMicrostructure phases => config_material%get('phase') allocate(counterPhase(phases%length),source=0) - homogenization => config_material%get('homogenization') - allocate(counterHomogenization(homogenization%length),source=0) + homogenizations => config_material%get('homogenization') + allocate(counterHomogenization(homogenizations%length),source=0) do e = 1, discretization_nElem microstructure => microstructures%get(discretization_microstructureAt(e)) From 79d672f4a77f53b9f1cb5d632daa7efc3a567a2e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Sep 2020 17:04:06 +0200 Subject: [PATCH 10/15] cleaning --- src/grid/DAMASK_grid.f90 | 4 ++-- src/grid/grid_mech_FEM.f90 | 1 - 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index f0589c3d7..ac30d487d 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -211,7 +211,7 @@ program DAMASK_grid if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo newLoadCase%deformation%mask = transpose(reshape(temp_maskVector,[ 3,3])) ! mask in 3x3 notation - newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation + newLoadCase%deformation%values = math_9to33(temp_valueVector) ! values in 3x3 notation case('p','stress', 's') temp_valueVector = 0.0_pReal do j = 1, 9 @@ -219,7 +219,7 @@ program DAMASK_grid if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) ! read value where applicable enddo newLoadCase%stress%mask = transpose(reshape(temp_maskVector,[ 3,3])) - newLoadCase%stress%values = math_9to33(temp_valueVector) + newLoadCase%stress%values = math_9to33(temp_valueVector) case('t','time','delta') ! increment time newLoadCase%time = IO_floatValue(line,chunkPos,i+1) case('n','incs','increments') ! number of increments diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 0fd7315a4..b680dafa9 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -62,7 +62,6 @@ module grid_mech_FEM real(pReal), dimension(3,3) :: & F_aimDot = 0.0_pReal, & !< assumed rate of average deformation gradient F_aim = math_I3, & !< current prescribed deformation gradient - F_aim_lastIter = math_I3, & F_aim_lastInc = math_I3, & !< previous average deformation gradient P_av = 0.0_pReal, & !< average 1st Piola--Kirchhoff stress P_aim = 0.0_pReal From 421d4b8f37c33103bc9c696879ac8c0340204f4e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Sep 2020 15:09:24 +0200 Subject: [PATCH 11/15] forgotten renames --- src/grid/DAMASK_grid.f90 | 2 +- src/material.f90 | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index ac30d487d..d092fb477 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -288,7 +288,7 @@ program DAMASK_grid errorID = 838 ! no rotation is allowed by stress BC print*, ' stress / GPa:' do i = 1, 3; do j = 1, 3 - if(newLoadCase%stress%maskLogical(i,j)) then + if(newLoadCase%stress%mask(i,j)) then write(IO_STDOUT,'(2x,f12.7)',advance='no') newLoadCase%stress%values(i,j)*1e-9_pReal else write(IO_STDOUT,'(2x,12a)',advance='no') ' * ' diff --git a/src/material.f90 b/src/material.f90 index 30c5f09a0..b6d043183 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -369,7 +369,7 @@ subroutine material_parseMicrostructure microstructure => microstructures%get(discretization_microstructureAt(e)) constituents => microstructure%get('constituents') - material_homogenizationAt(e) = homogenization%getIndex(microstructure%get_asString('homogenization')) + material_homogenizationAt(e) = homogenizations%getIndex(microstructure%get_asString('homogenization')) do i = 1, discretization_nIP counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) From 8396af5aec875916043edcd8510d600cf5d052ea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Sep 2020 15:19:31 +0200 Subject: [PATCH 12/15] not needed --- src/grid/DAMASK_grid.f90 | 4 ++-- src/grid/grid_damage_spectral.f90 | 6 ++---- src/grid/grid_thermal_spectral.f90 | 6 ++---- 3 files changed, 6 insertions(+), 10 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index d092fb477..df7df4d3b 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -431,9 +431,9 @@ program DAMASK_grid case(FIELD_MECH_ID) solres(field) = mech_solution(incInfo) case(FIELD_THERMAL_ID) - solres(field) = grid_thermal_spectral_solution(timeinc,timeIncOld) + solres(field) = grid_thermal_spectral_solution(timeinc) case(FIELD_DAMAGE_ID) - solres(field) = grid_damage_spectral_solution(timeinc,timeIncOld) + solres(field) = grid_damage_spectral_solution(timeinc) end select if (.not. solres(field)%converged) exit ! no solution found diff --git a/src/grid/grid_damage_spectral.f90 b/src/grid/grid_damage_spectral.f90 index 27dc5c1bd..7e529fc68 100644 --- a/src/grid/grid_damage_spectral.f90 +++ b/src/grid/grid_damage_spectral.f90 @@ -156,11 +156,10 @@ end subroutine grid_damage_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral damage scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) +function grid_damage_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment + timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull @@ -174,7 +173,6 @@ function grid_damage_spectral_solution(timeinc,timeinc_old) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc - params%timeincOld = timeinc_old call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 49be5ad7e..b4f9acddb 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -148,11 +148,10 @@ end subroutine grid_thermal_spectral_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the spectral thermal scheme with internal iterations !-------------------------------------------------------------------------------------------------- -function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) +function grid_thermal_spectral_solution(timeinc) result(solution) real(pReal), intent(in) :: & - timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment + timeinc !< increment in time for current solution integer :: i, j, k, cell type(tSolutionState) :: solution PetscInt :: devNull @@ -166,7 +165,6 @@ function grid_thermal_spectral_solution(timeinc,timeinc_old) result(solution) !-------------------------------------------------------------------------------------------------- ! set module wide availabe data params%timeinc = timeinc - params%timeincOld = timeinc_old call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) From 32b81770d9dcf043ad4dd4e9b1a416987c676ef5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Sep 2020 21:35:47 +0200 Subject: [PATCH 13/15] shorter --- src/grid/DAMASK_grid.f90 | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/grid/DAMASK_grid.f90 b/src/grid/DAMASK_grid.f90 index df7df4d3b..a1bb23375 100644 --- a/src/grid/DAMASK_grid.f90 +++ b/src/grid/DAMASK_grid.f90 @@ -365,11 +365,7 @@ program DAMASK_grid timeinc = loadCases(currentLoadCase)%time/real(loadCases(currentLoadCase)%incs,pReal) else if (currentLoadCase == 1) then ! 1st load case of logarithmic scale - if (inc == 1) then ! 1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real( 1-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st load case of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1-loadCases(1)%incs ,pReal)) - endif + timeinc = loadCases(1)%time*(2.0_pReal**real(max(inc-1,1)-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd else ! not-1st load case of logarithmic scale timeinc = time0 * & ( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc ,pReal)/& From 002fe04d353b94d70a2aecd832f28ebf77a31aef Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 26 Sep 2020 12:09:45 +0200 Subject: [PATCH 14/15] test for stress ramp --- PRIVATE | 2 +- src/grid/grid_mech_spectral_basic.f90 | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/PRIVATE b/PRIVATE index dc568df60..96a564660 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit dc568df60e36b659d9a1f84ac93fd4abb1b8fe3c +Subproject commit 96a564660b1f109009039785ca6cc4c670867178 diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index 8072f49b9..e7e544886 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -123,13 +123,13 @@ subroutine grid_mech_spectral_basic_init ! read numerical parameters and do sanity checks num_grid => config_numerics%get('grid',defaultVal=emptyDict) - num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) - num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) - num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) - num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) - num%eps_stress_rtol = num_grid%get_asFloat('eps_stress_rtol',defaultVal=1.0e-3_pReal) - num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) - num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) + num%update_gamma = num_grid%get_asBool ('update_gamma', defaultVal=.false.) + num%eps_div_atol = num_grid%get_asFloat ('eps_div_atol', defaultVal=1.0e-4_pReal) + num%eps_div_rtol = num_grid%get_asFloat ('eps_div_rtol', defaultVal=5.0e-4_pReal) + num%eps_stress_atol = num_grid%get_asFloat ('eps_stress_atol',defaultVal=1.0e3_pReal) + num%eps_stress_rtol = num_grid%get_asFloat ('eps_stress_rtol',defaultVal=1.0e-3_pReal) + num%itmin = num_grid%get_asInt ('itmin',defaultVal=1) + num%itmax = num_grid%get_asInt ('itmax',defaultVal=250) if (num%eps_div_atol <= 0.0_pReal) call IO_error(301,ext_msg='eps_div_atol') if (num%eps_div_rtol < 0.0_pReal) call IO_error(301,ext_msg='eps_div_rtol') From 4721743a57ac5b00a8cec00eec252fde8415eb40 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 26 Sep 2020 13:17:17 +0200 Subject: [PATCH 15/15] more thorough testing --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 96a564660..b1cb4c7f3 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 96a564660b1f109009039785ca6cc4c670867178 +Subproject commit b1cb4c7f306b3412704615793d6f61f4218ca24a