From a2e942033607a9d764aa5d64154001b364d7b5ff Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 14 Sep 2020 14:58:44 +0200 Subject: [PATCH 01/39] 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/39] 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/39] 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/39] 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/39] 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/39] 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/39] 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/39] 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/39] 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 c8dcd4b4ae3658877e3c37781b2578ff3ae1e5da Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 22 Sep 2020 23:51:47 +0200 Subject: [PATCH 10/39] 'decide' function should handle everything bug fixed in case of only flow yaml --- src/YAML_parse.f90 | 38 +++++++++++++++++++++++++------------- 1 file changed, 25 insertions(+), 13 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index d655bc2dc..6cefc5443 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -552,19 +552,17 @@ function to_flow(blck) s_flow, & !< start position in flow offset, & !< counts leading '- ' in nested lists end_line - if(isFlow(blck)) then - to_flow = trim(adjustl(blck)) - else - allocate(character(len=len(blck)*2)::to_flow) - ! move forward here (skip empty lines) and remove '----' if found - s_flow = 1 - s_blck = 1 - offset = 0 - call decide(blck,to_flow,s_blck,s_flow,offset) - to_flow = trim(to_flow(:s_flow-1)) - endif - end_line = index(to_flow,IO_EOL) - if(end_line > 0) to_flow = to_flow(:end_line-1) + + allocate(character(len=len(blck)*2)::to_flow) + ! move forward here (skip empty lines) and remove '----' if found + s_flow = 1 + s_blck = 1 + offset = 0 + call decide(blck,to_flow,s_blck,s_flow,offset) + to_flow = trim(to_flow(:s_flow-1)) + + end_line = index(to_flow,IO_EOL) + if(end_line > 0) to_flow = to_flow(:end_line-1) end function to_flow @@ -636,6 +634,20 @@ subroutine selfTest if (.not. to_flow(block_dict_newline) == flow_dict) error stop 'to_flow' end block basic_dict + only_flow: block + character(len=*), parameter :: flow_dict = & + " {a: [b,c: {d: e}, f: g, e]}"//IO_EOL + character(len=*), parameter :: flow_list = & + " [a,b: c, d,e: {f: g}]"//IO_EOL + character(len=*), parameter :: flow_1 = & + "{a: [b, {c: {d: e}}, {f: g}, e]}" + character(len=*), parameter :: flow_2 = & + "[a, {b: c}, d, {e: {f: g}}]" + + if (.not. to_flow(flow_dict) == flow_1) error stop 'to_flow' + if (.not. to_flow(flow_list) == flow_2) error stop 'to_flow' + end block only_flow + basic_flow: block character(len=*), parameter :: flow_braces = & " source: [{param: 1}, {param: 2}, {param: 3}, {param: 4}]"//IO_EOL From 79d672f4a77f53b9f1cb5d632daa7efc3a567a2e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 24 Sep 2020 17:04:06 +0200 Subject: [PATCH 11/39] 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 0de54404ee38f53cca8c6cd20921bf0641f059a8 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 25 Sep 2020 04:07:40 +0200 Subject: [PATCH 12/39] skip empty lines; yaml file optional start/stop indicator can be added --- src/IO.f90 | 2 ++ src/YAML_parse.f90 | 71 ++++++++++++++++++++++++++++++++++++++++------ 2 files changed, 64 insertions(+), 9 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index f9e708ead..eeb4ec8c6 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -494,6 +494,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'Unsupported feature' case (706) msg = 'Access by incorrect node type' + case (707) + msg = 'Abrupt end of file' !------------------------------------------------------------------------------------------------- ! errors related to the grid solver diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 6cefc5443..ade104cb2 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -227,6 +227,35 @@ logical function isKey(line) end function isKey +!-------------------------------------------------------------------------------------------------- +! @brief skip empty lines +! @details update start position in the block by skipping empty lines if present. +!-------------------------------------------------------------------------------------------------- +subroutine skip_empty_lines(blck,s_blck) + + character(len=*), intent(in) :: blck + integer, intent(inout) :: s_blck + + character(len=:), allocatable :: line + integer :: e_blck + logical :: not_empty + + not_empty = .false. + do while(.not. not_empty) + e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 + line = IO_rmComment(blck(s_blck:e_blck)) + if(len_trim(line) == 0) then + s_blck = e_blck + 2 + not_empty = .false. + else + not_empty = .true. + endif + enddo + + +end subroutine skip_empty_lines + + !-------------------------------------------------------------------------------------------------- ! @brief reads a line of YAML block which is already in flow style ! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser @@ -363,7 +392,9 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset) do while (s_blck <= len_trim(blck)) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if (len_trim(line) == 0) then + if(trim(line) == '---') then + exit + elseif (len_trim(line) == 0) then s_blck = e_blck + 2 ! forward to next line cycle elseif(indentDepth(line,offset) > indent) then @@ -377,8 +408,10 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset) else if(trim(adjustl(line)) == '-') then ! list item in next line s_blck = e_blck + 2 - e_blck = e_blck + index(blck(e_blck+2:),IO_EOL) + call skip_empty_lines(blck,s_blck) + e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) + if(trim(line) == '---') call IO_error(707,ext_msg=line) if(indentDepth(line) < indent .or. indentDepth(line) == indent) & call IO_error(701,ext_msg=line) @@ -447,7 +480,9 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset) do while (s_blck <= len_trim(blck)) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if (len_trim(line) == 0) then + if(trim(line) == '---') then + exit + elseif (len_trim(line) == 0) then s_blck = e_blck + 2 ! forward to next line cycle elseif(indentDepth(line,offset) < indent) then @@ -510,10 +545,12 @@ recursive subroutine decide(blck,flow,s_blck,s_flow,offset) character(len=:), allocatable :: line if(s_blck <= len(blck)) then + call skip_empty_lines(blck,s_blck) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - - if(len_trim(line) == 0) then + if(trim(line) == '---') then + continue ! end parsing at this point but not stop the simulation + elseif(len_trim(line) == 0) then s_blck = e_blck +2 call decide(blck,flow,s_blck,s_flow,offset) elseif (isListItem(line)) then @@ -548,16 +585,24 @@ function to_flow(blck) character(len=:), allocatable :: to_flow character(len=*), intent(in) :: blck !< YAML mixed style + + character(len=:), allocatable :: line integer :: s_blck, & !< start position in blck + e_blck, & !< end position of a line s_flow, & !< start position in flow offset, & !< counts leading '- ' in nested lists end_line - + allocate(character(len=len(blck)*2)::to_flow) - ! move forward here (skip empty lines) and remove '----' if found s_flow = 1 s_blck = 1 offset = 0 + + call skip_empty_lines(blck,s_blck) + e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 + line = IO_rmComment(blck(s_blck:e_blck)) + if(trim(line) == '---') s_blck = e_blck + 2 + call decide(blck,to_flow,s_blck,s_flow,offset) to_flow = trim(to_flow(:s_flow-1)) @@ -662,12 +707,20 @@ subroutine selfTest basic_mixed: block character(len=*), parameter :: block_flow = & + " "//IO_EOL//& + " "//IO_EOL//& + "---"//IO_EOL//& " aa:"//IO_EOL//& " - "//IO_EOL//& - " param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//& + " "//IO_EOL//& + " "//IO_EOL//& + " param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//& " - c: d"//IO_EOL//& " bb:"//IO_EOL//& - " - {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL + " "//IO_EOL//& + " - "//IO_EOL//& + " {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL//& + "---"//IO_EOL character(len=*), parameter :: mixed_flow = & "{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}" From 21ff587e1720c53eb7bf441860e05a29f3f7ed8d Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 25 Sep 2020 04:52:03 +0200 Subject: [PATCH 13/39] better logic --- src/YAML_parse.f90 | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index ade104cb2..218d4688b 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -238,18 +238,14 @@ subroutine skip_empty_lines(blck,s_blck) character(len=:), allocatable :: line integer :: e_blck - logical :: not_empty + logical :: empty - not_empty = .false. - do while(.not. not_empty) + empty = .true. + do while(empty) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if(len_trim(line) == 0) then - s_blck = e_blck + 2 - not_empty = .false. - else - not_empty = .true. - endif + empty = len_trim(line) == 0 + if(empty) s_blck = e_blck + 2 enddo From 08f5851c82771fb4365d898c62a104558ddf6d29 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Fri, 25 Sep 2020 10:59:03 +0200 Subject: [PATCH 14/39] take care of empty lines in this slightly new setup --- src/YAML_parse.f90 | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 218d4688b..65b9afe3d 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -594,12 +594,13 @@ function to_flow(blck) s_blck = 1 offset = 0 - call skip_empty_lines(blck,s_blck) - e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 - line = IO_rmComment(blck(s_blck:e_blck)) - if(trim(line) == '---') s_blck = e_blck + 2 - - call decide(blck,to_flow,s_blck,s_flow,offset) + if(len_trim(blck) /= 0) then + call skip_empty_lines(blck,s_blck) + e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 + line = IO_rmComment(blck(s_blck:e_blck)) + if(trim(line) == '---') s_blck = e_blck + 2 + call decide(blck,to_flow,s_blck,s_flow,offset) + endif to_flow = trim(to_flow(:s_flow-1)) end_line = index(to_flow,IO_EOL) From 421d4b8f37c33103bc9c696879ac8c0340204f4e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 25 Sep 2020 15:09:24 +0200 Subject: [PATCH 15/39] 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 16/39] 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 17/39] 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 18/39] 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 19/39] 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 From c46b4d90a613a2ad881d9eec578326fc21e918ed Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 12:48:29 +0200 Subject: [PATCH 20/39] modularizing --- src/crystallite.f90 | 22 +++++++++++++++++----- 1 file changed, 17 insertions(+), 5 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 0bd816142..c0c140571 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1241,7 +1241,7 @@ subroutine integrateStateFPI(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck contains @@ -1328,7 +1328,7 @@ subroutine integrateStateEuler(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck end subroutine integrateStateEuler @@ -1422,7 +1422,7 @@ subroutine integrateStateAdaptiveEuler(todo) enddo; enddo; enddo !$OMP END PARALLEL DO - if (nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck end subroutine integrateStateAdaptiveEuler @@ -1612,7 +1612,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB) enddo; enddo; enddo !$OMP END PARALLEL DO - if(nonlocalBroken) call nonlocalConvergenceCheck + call nonlocalConvergenceCheck end subroutine integrateStateRK @@ -1624,7 +1624,19 @@ end subroutine integrateStateRK subroutine nonlocalConvergenceCheck integer :: e,i,p - + logical :: nonlocal_broken + + nonlocal_broken = .false. + !$OMP PARALLEL DO PRIVATE(p) + do e = FEsolving_execElem(1),FEsolving_execElem(2) + p = material_phaseAt(1,e) + do i = FEsolving_execIP(1),FEsolving_execIP(2) + if(plasticState(p)%nonlocal .and. .not. crystallite_converged(1,i,e)) nonlocal_broken = .true. + enddo + enddo + !$OMP END PARALLEL DO + + if(.not. nonlocal_broken) return !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) p = material_phaseAt(1,e) From f1e96489cc3d59fe11420a9d395caf025fb234f1 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 12:56:12 +0200 Subject: [PATCH 21/39] better readable somehow on the cost of the nonlocal performance --- src/crystallite.f90 | 43 ++++++++++++++++--------------------------- 1 file changed, 16 insertions(+), 27 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index c0c140571..5f4b3bbf7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1147,16 +1147,15 @@ subroutine integrateStateFPI(todo) plastic_dotState real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState logical :: & - nonlocalBroken, broken + broken + - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then - + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1164,7 +1163,6 @@ subroutine integrateStateFPI(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle size_pl = plasticState(p)%sizeDotState @@ -1236,7 +1234,6 @@ subroutine integrateStateFPI(todo) endif enddo iteration - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. endif enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1284,16 +1281,15 @@ subroutine integrateStateEuler(todo) s, & sizeDotState logical :: & - nonlocalBroken, broken + broken - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1301,7 +1297,6 @@ subroutine integrateStateEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle sizeDotState = plasticState(p)%sizeDotState @@ -1318,11 +1313,9 @@ subroutine integrateStateEuler(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle broken = integrateStress(g,i,e) - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. crystallite_converged(g,i,e) = .not. broken endif enddo; enddo; enddo @@ -1349,20 +1342,19 @@ subroutine integrateStateAdaptiveEuler(todo) s, & sizeDotState logical :: & - nonlocalBroken, broken + broken real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) broken = .false. - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1418,7 +1410,6 @@ subroutine integrateStateAdaptiveEuler(todo) enddo endif - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1503,19 +1494,18 @@ subroutine integrateStateRK(todo,A,B,CC,DB) s, & sizeDotState logical :: & - nonlocalBroken, broken + broken real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - nonlocalBroken = .false. !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken) do e = FEsolving_execElem(1),FEsolving_execElem(2) do i = FEsolving_execIP(1),FEsolving_execIP(2) do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) broken = .false. - p = material_phaseAt(g,e) - if(todo(g,i,e) .and. .not. (nonlocalBroken .and. plasticState(p)%nonlocal)) then + if(todo(g,i,e)) then + p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & @@ -1608,7 +1598,6 @@ subroutine integrateStateRK(todo,A,B,CC,DB) crystallite_converged(g,i,e) = .not. broken endif - if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. enddo; enddo; enddo !$OMP END PARALLEL DO @@ -1624,8 +1613,8 @@ end subroutine integrateStateRK subroutine nonlocalConvergenceCheck integer :: e,i,p - logical :: nonlocal_broken - + logical :: nonlocal_broken + nonlocal_broken = .false. !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) @@ -1635,7 +1624,7 @@ subroutine nonlocalConvergenceCheck enddo enddo !$OMP END PARALLEL DO - + if(.not. nonlocal_broken) return !$OMP PARALLEL DO PRIVATE(p) do e = FEsolving_execElem(1),FEsolving_execElem(2) From 587d5ee445394808eb55ad2e9164575428e589ea Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 13:13:53 +0200 Subject: [PATCH 22/39] no need for two loops --- src/crystallite.f90 | 97 ++++++++++++--------------------------------- 1 file changed, 25 insertions(+), 72 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 5f4b3bbf7..d21fd453f 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -460,16 +460,18 @@ function crystallite_stress() math_inv33(crystallite_Fi(1:3,1:3,c,i,e))) crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) crystallite_converged(c,i,e) = .false. + call integrateState(c,i,e) endif enddo enddo enddo elementLooping3 !$OMP END PARALLEL DO + + call nonlocalConvergenceCheck !-------------------------------------------------------------------------------------------------- ! integrate --- requires fully defined state array (basic + dependent state) - if (any(todo)) call integrateState(todo) ! TODO: unroll into proper elementloop to avoid N^2 for single point evaluation where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation @@ -1125,9 +1127,8 @@ end function integrateStress !> @brief integrate stress, state with adaptive 1st order explicit Euler method !> using Fixed Point Iteration to adapt the stepsize !-------------------------------------------------------------------------------------------------- -subroutine integrateStateFPI(todo) +subroutine integrateStateFPI(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo integer :: & NiterationState, & !< number of iterations in state loop e, & !< element index in element loop @@ -1150,11 +1151,6 @@ subroutine integrateStateFPI(todo) broken - !$OMP PARALLEL DO PRIVATE(size_pl,size_so,r,zeta,p,c,plastic_dotState,source_dotState,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1163,7 +1159,7 @@ subroutine integrateStateFPI(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return size_pl = plasticState(p)%sizeDotState plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & @@ -1234,11 +1230,7 @@ subroutine integrateStateFPI(todo) endif enddo iteration - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - call nonlocalConvergenceCheck contains @@ -1268,9 +1260,7 @@ end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- !> @brief integrate state with 1st order explicit Euler method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateEuler(todo) - - logical, dimension(:,:,:), intent(in) :: todo +subroutine integrateStateEuler(g,i,e) integer :: & e, & !< element index in element loop @@ -1283,11 +1273,6 @@ subroutine integrateStateEuler(todo) logical :: & broken - !$OMP PARALLEL DO PRIVATE (sizeDotState,p,c,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1297,7 +1282,7 @@ subroutine integrateStateEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & @@ -1313,15 +1298,10 @@ subroutine integrateStateEuler(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + if(broken) return broken = integrateStress(g,i,e) crystallite_converged(g,i,e) = .not. broken - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - call nonlocalConvergenceCheck end subroutine integrateStateEuler @@ -1329,9 +1309,7 @@ end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- !> @brief integrate stress, state with 1st order Euler method with adaptive step size !-------------------------------------------------------------------------------------------------- -subroutine integrateStateAdaptiveEuler(todo) - - logical, dimension(:,:,:), intent(in) :: todo +subroutine integrateStateAdaptiveEuler(g,i,e) integer :: & e, & ! element index in element loop @@ -1347,13 +1325,7 @@ subroutine integrateStateAdaptiveEuler(todo) real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,residuum_plastic,residuum_source,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - broken = .false. - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1362,7 +1334,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState @@ -1381,17 +1353,17 @@ subroutine integrateStateAdaptiveEuler(todo) broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + if(broken) return broken = integrateStress(g,i,e) - if(broken) cycle + if(broken) return broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState @@ -1407,13 +1379,7 @@ subroutine integrateStateAdaptiveEuler(todo) + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - call nonlocalConvergenceCheck + enddo end subroutine integrateStateAdaptiveEuler @@ -1421,9 +1387,9 @@ end subroutine integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the classic Runge Kutta method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRK4(todo) +subroutine integrateStateRK4(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo + integer :: g,i,e real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1436,7 +1402,7 @@ subroutine integrateStateRK4(todo) real(pReal), dimension(4), parameter :: & B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal] - call integrateStateRK(todo,A,B,C) + call integrateStateRK(g,i,e,A,B,C) end subroutine integrateStateRK4 @@ -1444,9 +1410,9 @@ end subroutine integrateStateRK4 !--------------------------------------------------------------------------------------------------- !> @brief Integrate state (including stress integration) with the Cash-Carp method !--------------------------------------------------------------------------------------------------- -subroutine integrateStateRKCK45(todo) +subroutine integrateStateRKCK45(g,i,e) - logical, dimension(:,:,:), intent(in) :: todo + integer :: g,i,e real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1466,7 +1432,7 @@ subroutine integrateStateRKCK45(todo) [2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,& 13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal] - call integrateStateRK(todo,A,B,C,DB) + call integrateStateRK(g,i,e,A,B,C,DB) end subroutine integrateStateRKCK45 @@ -1475,9 +1441,8 @@ end subroutine integrateStateRKCK45 !> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an !! embedded explicit Runge-Kutta method !-------------------------------------------------------------------------------------------------- -subroutine integrateStateRK(todo,A,B,CC,DB) +subroutine integrateStateRK(g,i,e,A,B,CC,DB) - logical, dimension(:,:,:), intent(in) :: todo real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(:), intent(in) :: B, CC @@ -1498,13 +1463,6 @@ subroutine integrateStateRK(todo,A,B,CC,DB) real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - !$OMP PARALLEL DO PRIVATE(sizeDotState,p,c,plastic_RKdotState,source_RKdotState,broken) - do e = FEsolving_execElem(1),FEsolving_execElem(2) - do i = FEsolving_execIP(1),FEsolving_execIP(2) - do g = 1,homogenization_Ngrains(material_homogenizationAt(e)) - broken = .false. - - if(todo(g,i,e)) then p = material_phaseAt(g,e) c = material_phaseMemberAt(g,i,e) @@ -1513,7 +1471,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB) crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) cycle + if(broken) return do stage = 1,size(A,1) sizeDotState = plasticState(p)%sizeDotState @@ -1558,7 +1516,7 @@ subroutine integrateStateRK(todo,A,B,CC,DB) if(broken) exit enddo - if(broken) cycle + if(broken) return sizeDotState = plasticState(p)%sizeDotState @@ -1587,21 +1545,16 @@ subroutine integrateStateRK(todo,A,B,CC,DB) sourceState(p)%p(s)%state(1:sizeDotState,c), & sourceState(p)%p(s)%atol(1:sizeDotState)) enddo - if(broken) cycle + if(broken) return broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) cycle + if(broken) return broken = integrateStress(g,i,e) crystallite_converged(g,i,e) = .not. broken - endif - enddo; enddo; enddo - !$OMP END PARALLEL DO - - call nonlocalConvergenceCheck end subroutine integrateStateRK From d0df748fc142481f29cdda9232634aade01b56f8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 13:36:15 +0200 Subject: [PATCH 23/39] cleaning --- src/crystallite.f90 | 474 ++++++++++++++++++++++---------------------- 1 file changed, 238 insertions(+), 236 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index d21fd453f..9411cdfc7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -467,7 +467,7 @@ function crystallite_stress() enddo enddo elementLooping3 !$OMP END PARALLEL DO - + call nonlocalConvergenceCheck !-------------------------------------------------------------------------------------------------- @@ -1129,11 +1129,12 @@ end function integrateStress !-------------------------------------------------------------------------------------------------- subroutine integrateStateFPI(g,i,e) - integer :: & - NiterationState, & !< number of iterations in state loop + integer, intent(in) :: & e, & !< element index in element loop i, & !< integration point index in ip loop - g, & !< grain index in grain loop + g !< grain index in grain loop + integer :: & + NiterationState, & !< number of iterations in state loop p, & c, & s, & @@ -1150,86 +1151,85 @@ subroutine integrateStateFPI(g,i,e) logical :: & broken + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + size_pl = plasticState(p)%sizeDotState + plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & + + plasticState(p)%dotState (1:size_pl,c) & + * crystallite_subdt(g,i,e) + plastic_dotState(1:size_pl,2) = 0.0_pReal + do s = 1, phase_Nsources(p) + size_so(s) = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & + + sourceState(p)%p(s)%dotState (1:size_so(s),c) & + * crystallite_subdt(g,i,e) + source_dotState(1:size_so(s),2,s) = 0.0_pReal + enddo - size_pl = plasticState(p)%sizeDotState - plasticState(p)%state(1:size_pl,c) = plasticState(p)%subState0(1:size_pl,c) & - + plasticState(p)%dotState (1:size_pl,c) & - * crystallite_subdt(g,i,e) - plastic_dotState(1:size_pl,2) = 0.0_pReal - do s = 1, phase_Nsources(p) - size_so(s) = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%subState0(1:size_so(s),c) & - + sourceState(p)%p(s)%dotState (1:size_so(s),c) & - * crystallite_subdt(g,i,e) - source_dotState(1:size_so(s),2,s) = 0.0_pReal - enddo + iteration: do NiterationState = 1, num%nState - iteration: do NiterationState = 1, num%nState + if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) + plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) + do s = 1, phase_Nsources(p) + if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) + source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) + enddo - if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1) - plastic_dotState(1:size_pl,1) = plasticState(p)%dotState(:,c) - do s = 1, phase_Nsources(p) - if(nIterationState > 1) source_dotState(1:size_so(s),2,s) = source_dotState(1:size_so(s),1,s) - source_dotState(1:size_so(s),1,s) = sourceState(p)%p(s)%dotState(:,c) - enddo + broken = integrateStress(g,i,e) + if(broken) exit iteration - broken = integrateStress(g,i,e) - if(broken) exit iteration + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) exit iteration - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) exit iteration + zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& + plastic_dotState(1:size_pl,2)) + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & + + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) + r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & + - plasticState(p)%subState0(1:size_pl,c) & + - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) + plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & + - r(1:size_pl) + crystallite_converged(g,i,e) = converged(r(1:size_pl), & + plasticState(p)%state(1:size_pl,c), & + plasticState(p)%atol(1:size_pl)) + do s = 1, phase_Nsources(p) + zeta = damper(sourceState(p)%p(s)%dotState(:,c), & + source_dotState(1:size_so(s),1,s),& + source_dotState(1:size_so(s),2,s)) + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & + + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) + r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & + - sourceState(p)%p(s)%subState0(1:size_so(s),c) & + - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & + - r(1:size_so(s)) + crystallite_converged(g,i,e) = & + crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & + sourceState(p)%p(s)%state(1:size_so(s),c), & + sourceState(p)%p(s)%atol(1:size_so(s))) + enddo - zeta = damper(plasticState(p)%dotState(:,c),plastic_dotState(1:size_pl,1),& - plastic_dotState(1:size_pl,2)) - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) * zeta & - + plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta) - r(1:size_pl) = plasticState(p)%state (1:size_pl,c) & - - plasticState(p)%subState0(1:size_pl,c) & - - plasticState(p)%dotState (1:size_pl,c) * crystallite_subdt(g,i,e) - plasticState(p)%state(1:size_pl,c) = plasticState(p)%state(1:size_pl,c) & - - r(1:size_pl) - crystallite_converged(g,i,e) = converged(r(1:size_pl), & - plasticState(p)%state(1:size_pl,c), & - plasticState(p)%atol(1:size_pl)) - do s = 1, phase_Nsources(p) - zeta = damper(sourceState(p)%p(s)%dotState(:,c), & - source_dotState(1:size_so(s),1,s),& - source_dotState(1:size_so(s),2,s)) - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) * zeta & - + source_dotState(1:size_so(s),1,s)* (1.0_pReal - zeta) - r(1:size_so(s)) = sourceState(p)%p(s)%state (1:size_so(s),c) & - - sourceState(p)%p(s)%subState0(1:size_so(s),c) & - - sourceState(p)%p(s)%dotState (1:size_so(s),c) * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:size_so(s),c) = sourceState(p)%p(s)%state(1:size_so(s),c) & - - r(1:size_so(s)) - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(r(1:size_so(s)), & - sourceState(p)%p(s)%state(1:size_so(s),c), & - sourceState(p)%p(s)%atol(1:size_so(s))) - enddo + if(crystallite_converged(g,i,e)) then + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + exit iteration + endif - if(crystallite_converged(g,i,e)) then - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - exit iteration - endif - - enddo iteration + enddo iteration contains @@ -1262,10 +1262,11 @@ end subroutine integrateStateFPI !-------------------------------------------------------------------------------------------------- subroutine integrateStateEuler(g,i,e) - integer :: & + integer, intent(in) :: & e, & !< element index in element loop i, & !< integration point index in ip loop - g, & !< grain index in grain loop + g !< grain index in grain loop + integer :: & p, & c, & s, & @@ -1273,35 +1274,34 @@ subroutine integrateStateEuler(g,i,e) logical :: & broken + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + enddo - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return - - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken end subroutine integrateStateEuler @@ -1311,10 +1311,11 @@ end subroutine integrateStateEuler !-------------------------------------------------------------------------------------------------- subroutine integrateStateAdaptiveEuler(g,i,e) + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop integer :: & - e, & ! element index in element loop - i, & ! integration point index in ip loop - g, & ! grain index in grain loop p, & c, & s, & @@ -1326,60 +1327,60 @@ subroutine integrateStateAdaptiveEuler(g,i,e) real(pReal), dimension(constitutive_source_maxSizeDotState,maxval(phase_Nsources)) :: residuum_source - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - sizeDotState = plasticState(p)%sizeDotState + sizeDotState = plasticState(p)%sizeDotState - residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState + residuum_plastic(1:sizeDotState) = - plasticState(p)%dotstate(1:sizeDotState,c) * 0.5_pReal * crystallite_subdt(g,i,e) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState - residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & - * 0.5_pReal * crystallite_subdt(g,i,e) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) - enddo + residuum_source(1:sizeDotState,s) = - sourceState(p)%p(s)%dotstate(1:sizeDotState,c) & + * 0.5_pReal * crystallite_subdt(g,i,e) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) + enddo - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return - broken = integrateStress(g,i,e) - if(broken) return + broken = integrateStress(g,i,e) + if(broken) return - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - sizeDotState = plasticState(p)%sizeDotState - crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & - + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + sizeDotState = plasticState(p)%sizeDotState + crystallite_converged(g,i,e) = converged(residuum_plastic(1:sizeDotState) & + + 0.5_pReal * plasticState(p)%dotState(:,c) * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - crystallite_converged(g,i,e) = & - crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & - + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + crystallite_converged(g,i,e) = & + crystallite_converged(g,i,e) .and. converged(residuum_source(1:sizeDotState,s) & + + 0.5_pReal*sourceState(p)%p(s)%dotState(:,c)*crystallite_subdt(g,i,e), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & + sourceState(p)%p(s)%atol(1:sizeDotState)) + enddo end subroutine integrateStateAdaptiveEuler @@ -1389,7 +1390,7 @@ end subroutine integrateStateAdaptiveEuler !--------------------------------------------------------------------------------------------------- subroutine integrateStateRK4(g,i,e) - integer :: g,i,e + integer, intent(in) :: g,i,e real(pReal), dimension(3,3), parameter :: & A = reshape([& @@ -1412,7 +1413,7 @@ end subroutine integrateStateRK4 !--------------------------------------------------------------------------------------------------- subroutine integrateStateRKCK45(g,i,e) - integer :: g,i,e + integer, intent(in) :: g,i,e real(pReal), dimension(5,5), parameter :: & A = reshape([& @@ -1448,10 +1449,11 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB) real(pReal), dimension(:), intent(in) :: B, CC real(pReal), dimension(:), intent(in), optional :: DB + integer, intent(in) :: & + e, & !< element index in element loop + i, & !< integration point index in ip loop + g !< grain index in grain loop integer :: & - e, & ! element index in element loop - i, & ! integration point index in ip loop - g, & ! grain index in grain loop stage, & ! stage index in integration stage loop n, & p, & @@ -1463,97 +1465,97 @@ subroutine integrateStateRK(g,i,e,A,B,CC,DB) real(pReal), dimension(constitutive_source_maxSizeDotState,size(B),maxval(phase_Nsources)) :: source_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState - p = material_phaseAt(g,e) - c = material_phaseMemberAt(g,i,e) + p = material_phaseAt(g,e) + c = material_phaseMemberAt(g,i,e) - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e,p,c) - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e), g,i,e,p,c) + if(broken) return - do stage = 1,size(A,1) - sizeDotState = plasticState(p)%sizeDotState - plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) - plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) - enddo + do stage = 1,size(A,1) + sizeDotState = plasticState(p)%sizeDotState + plastic_RKdotState(1:sizeDotState,stage) = plasticState(p)%dotState(:,c) + plasticState(p)%dotState(:,c) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + source_RKdotState(1:sizeDotState,stage,s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = A(1,stage) * source_RKdotState(1:sizeDotState,1,s) + enddo - do n = 2, stage - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & - + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & - + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) - enddo - enddo + do n = 2, stage + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%dotState(:,c) = plasticState(p)%dotState(:,c) & + + A(n,stage) * plastic_RKdotState(1:sizeDotState,n) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%dotState(:,c) = sourceState(p)%p(s)%dotState(:,c) & + + A(n,stage) * source_RKdotState(1:sizeDotState,n,s) + enddo + enddo - sizeDotState = plasticState(p)%sizeDotState - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - enddo - - broken = integrateStress(g,i,e,CC(stage)) - if(broken) exit - - broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_partionedF0, & - crystallite_Fi(1:3,1:3,g,i,e), & - crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) - if(broken) exit - - enddo - if(broken) return - - sizeDotState = plasticState(p)%sizeDotState - - plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) - plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) - plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & - + plasticState(p)%dotState (1:sizeDotState,c) & + sizeDotState = plasticState(p)%sizeDotState + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & - * crystallite_subdt(g,i,e), & - plasticState(p)%state(1:sizeDotState,c), & - plasticState(p)%atol(1:sizeDotState)) + enddo - do s = 1, phase_Nsources(p) - sizeDotState = sourceState(p)%p(s)%sizeDotState + broken = integrateStress(g,i,e,CC(stage)) + if(broken) exit - source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) - sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) - sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & - + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & - * crystallite_subdt(g,i,e) - if(present(DB)) & - broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),s),DB) & - * crystallite_subdt(g,i,e), & - sourceState(p)%p(s)%state(1:sizeDotState,c), & - sourceState(p)%p(s)%atol(1:sizeDotState)) - enddo - if(broken) return + broken = constitutive_collectDotState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_partionedF0, & + crystallite_Fi(1:3,1:3,g,i,e), & + crystallite_partionedFp0, & + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) + if(broken) exit - broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & - crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) - if(broken) return + enddo + if(broken) return - broken = integrateStress(g,i,e) - crystallite_converged(g,i,e) = .not. broken + sizeDotState = plasticState(p)%sizeDotState + + plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (p)%dotState(:,c) + plasticState(p)%dotState(:,c) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B) + plasticState(p)%state(1:sizeDotState,c) = plasticState(p)%subState0(1:sizeDotState,c) & + + plasticState(p)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + if(present(DB)) & + broken = .not. converged( matmul(plastic_RKdotState(1:sizeDotState,1:size(DB)),DB) & + * crystallite_subdt(g,i,e), & + plasticState(p)%state(1:sizeDotState,c), & + plasticState(p)%atol(1:sizeDotState)) + + do s = 1, phase_Nsources(p) + sizeDotState = sourceState(p)%p(s)%sizeDotState + + source_RKdotState(1:sizeDotState,size(B),s) = sourceState(p)%p(s)%dotState(:,c) + sourceState(p)%p(s)%dotState(:,c) = matmul(source_RKdotState(1:sizeDotState,1:size(B),s),B) + sourceState(p)%p(s)%state(1:sizeDotState,c) = sourceState(p)%p(s)%subState0(1:sizeDotState,c) & + + sourceState(p)%p(s)%dotState (1:sizeDotState,c) & + * crystallite_subdt(g,i,e) + if(present(DB)) & + broken = broken .or. .not. converged(matmul(source_RKdotState(1:sizeDotState,1:size(DB),s),DB) & + * crystallite_subdt(g,i,e), & + sourceState(p)%p(s)%state(1:sizeDotState,c), & + sourceState(p)%p(s)%atol(1:sizeDotState)) + enddo + if(broken) return + + broken = constitutive_deltaState(crystallite_S(1:3,1:3,g,i,e), & + crystallite_Fe(1:3,1:3,g,i,e), & + crystallite_Fi(1:3,1:3,g,i,e),g,i,e,p,c) + if(broken) return + + broken = integrateStress(g,i,e) + crystallite_converged(g,i,e) = .not. broken end subroutine integrateStateRK From a61bf3bb26d4eac3a04002225ac0cda10e832727 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 29 Sep 2020 19:33:30 +0200 Subject: [PATCH 24/39] file endings, file header can be added, take care of EOF --- src/IO.f90 | 5 +++++ src/YAML_parse.f90 | 50 ++++++++++++++++++++++++++++++++++------------ 2 files changed, 42 insertions(+), 13 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index eeb4ec8c6..0542e7a62 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -496,6 +496,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'Access by incorrect node type' case (707) msg = 'Abrupt end of file' + case (708) + msg = '--- expected after YAML file header' !------------------------------------------------------------------------------------------------- ! errors related to the grid solver @@ -623,6 +625,9 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg) msg = 'polar decomposition failed' case (700) msg = 'unknown crystal symmetry' + case (709) + msg = 'read only the first document' + case (850) msg = 'max number of cut back exceeded, terminating' case default diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index 65b9afe3d..bb990fc52 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -236,15 +236,13 @@ subroutine skip_empty_lines(blck,s_blck) character(len=*), intent(in) :: blck integer, intent(inout) :: s_blck - character(len=:), allocatable :: line integer :: e_blck logical :: empty empty = .true. - do while(empty) + do while(empty .and. len_trim(blck(s_blck:)) /= 0) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 - line = IO_rmComment(blck(s_blck:e_blck)) - empty = len_trim(line) == 0 + empty = len_trim(IO_rmComment(blck(s_blck:e_blck))) == 0 if(empty) s_blck = e_blck + 2 enddo @@ -252,6 +250,31 @@ subroutine skip_empty_lines(blck,s_blck) end subroutine skip_empty_lines +!-------------------------------------------------------------------------------------------------- +! @brief skip file header +! @details update start position in the block by skipping file header if present. +!-------------------------------------------------------------------------------------------------- +subroutine skip_file_header(blck,s_blck) + + character(len=*), intent(in) :: blck + integer, intent(inout) :: s_blck + + character(len=:), allocatable :: line + + line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2)) + if(adjustl(line(1:5)) == '%YAML') then + s_blck = s_blck + index(blck(s_blck:),IO_EOL) + call skip_empty_lines(blck,s_blck) + if(trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == '---') then + s_blck = s_blck + index(blck(s_blck:),IO_EOL) + else + call IO_error(708,ext_msg = line) + endif + endif + +end subroutine skip_file_header + + !-------------------------------------------------------------------------------------------------- ! @brief reads a line of YAML block which is already in flow style ! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser @@ -388,7 +411,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset) do while (s_blck <= len_trim(blck)) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if(trim(line) == '---') then + if(trim(line) == '---' .or. trim(line) == '...') then exit elseif (len_trim(line) == 0) then s_blck = e_blck + 2 ! forward to next line @@ -476,7 +499,7 @@ recursive subroutine dct(blck,flow,s_blck,s_flow,offset) do while (s_blck <= len_trim(blck)) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if(trim(line) == '---') then + if(trim(line) == '---' .or. trim(line) == '...') then exit elseif (len_trim(line) == 0) then s_blck = e_blck + 2 ! forward to next line @@ -544,7 +567,7 @@ recursive subroutine decide(blck,flow,s_blck,s_flow,offset) call skip_empty_lines(blck,s_blck) e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 line = IO_rmComment(blck(s_blck:e_blck)) - if(trim(line) == '---') then + if(trim(line) == '---' .or. trim(line) == '...') then continue ! end parsing at this point but not stop the simulation elseif(len_trim(line) == 0) then s_blck = e_blck +2 @@ -584,7 +607,6 @@ function to_flow(blck) character(len=:), allocatable :: line integer :: s_blck, & !< start position in blck - e_blck, & !< end position of a line s_flow, & !< start position in flow offset, & !< counts leading '- ' in nested lists end_line @@ -596,13 +618,14 @@ function to_flow(blck) if(len_trim(blck) /= 0) then call skip_empty_lines(blck,s_blck) - e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 - line = IO_rmComment(blck(s_blck:e_blck)) - if(trim(line) == '---') s_blck = e_blck + 2 + call skip_file_header(blck,s_blck) + line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2)) + if(trim(line) == '---') s_blck = s_blck + index(blck(s_blck:),IO_EOL) call decide(blck,to_flow,s_blck,s_flow,offset) endif + line = IO_rmComment(blck(s_blck:s_blck+index(blck(s_blck:),IO_EOL)-2)) + if(trim(line)== '---') call IO_warning(709,ext_msg=line) to_flow = trim(to_flow(:s_flow-1)) - end_line = index(to_flow,IO_EOL) if(end_line > 0) to_flow = to_flow(:end_line-1) @@ -704,6 +727,7 @@ subroutine selfTest basic_mixed: block character(len=*), parameter :: block_flow = & + "%YAML 1.2"//IO_EOL//& " "//IO_EOL//& " "//IO_EOL//& "---"//IO_EOL//& @@ -717,7 +741,7 @@ subroutine selfTest " "//IO_EOL//& " - "//IO_EOL//& " {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL//& - "---"//IO_EOL + "..."//IO_EOL character(len=*), parameter :: mixed_flow = & "{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}" From 45b906906d90bfec363f63598c29c4dae0d3ed4e Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 29 Sep 2020 20:07:33 +0200 Subject: [PATCH 25/39] test before reading config files --- src/CPFEM.f90 | 4 ++-- src/CPFEM2.f90 | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index 05255e2a1..f93472e51 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -78,11 +78,11 @@ subroutine CPFEM_initAll call DAMASK_interface_init call prec_init call IO_init + call YAML_types_init + call YAML_parse_init call config_init call math_init call rotations_init - call YAML_types_init - call YAML_parse_init call HDF5_utilities_init call results_init(.false.) call discretization_marc_init diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index fed43ba78..994859758 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -48,11 +48,11 @@ subroutine CPFEM_initAll #ifdef Mesh call FEM_quadrature_init #endif + call YAML_types_init + call YAML_parse_init call config_init call math_init call rotations_init - call YAML_types_init - call YAML_parse_init call lattice_init call HDF5_utilities_init call results_init(restart=interface_restartInc>0) From 385cda9224dba4e84d0e93f5aba967efed279811 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 29 Sep 2020 20:13:02 +0200 Subject: [PATCH 26/39] remove unnecessary variables --- src/YAML_parse.f90 | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index bb990fc52..a63a07cf1 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -236,14 +236,12 @@ subroutine skip_empty_lines(blck,s_blck) character(len=*), intent(in) :: blck integer, intent(inout) :: s_blck - integer :: e_blck logical :: empty empty = .true. do while(empty .and. len_trim(blck(s_blck:)) /= 0) - e_blck = s_blck + index(blck(s_blck:),IO_EOL) - 2 - empty = len_trim(IO_rmComment(blck(s_blck:e_blck))) == 0 - if(empty) s_blck = e_blck + 2 + empty = len_trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == 0 + if(empty) s_blck = s_blck + index(blck(s_blck:),IO_EOL) enddo From 15ef6c8ceb65e671b6c950af4ce62c153da47421 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Wed, 30 Sep 2020 01:50:10 +0200 Subject: [PATCH 27/39] more fortran like --- src/YAML_parse.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/YAML_parse.f90 b/src/YAML_parse.f90 index a63a07cf1..cb5b726dc 100644 --- a/src/YAML_parse.f90 +++ b/src/YAML_parse.f90 @@ -260,7 +260,7 @@ subroutine skip_file_header(blck,s_blck) character(len=:), allocatable :: line line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2)) - if(adjustl(line(1:5)) == '%YAML') then + if(index(adjustl(line),'%YAML') == 1) then s_blck = s_blck + index(blck(s_blck:),IO_EOL) call skip_empty_lines(blck,s_blck) if(trim(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))) == '---') then @@ -725,7 +725,7 @@ subroutine selfTest basic_mixed: block character(len=*), parameter :: block_flow = & - "%YAML 1.2"//IO_EOL//& + "%YAML 1.1"//IO_EOL//& " "//IO_EOL//& " "//IO_EOL//& "---"//IO_EOL//& From bd5d557fbb560d08ace9d3fa8b0d3c91832d7639 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 29 Sep 2020 23:32:17 +0200 Subject: [PATCH 28/39] untangling the spaghetti subLp and subLi are local variables --- src/crystallite.f90 | 50 +++++++++++++-------------------------------- 1 file changed, 14 insertions(+), 36 deletions(-) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 9411cdfc7..ad1c2d4a7 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -54,12 +54,10 @@ module crystallite ! crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc crystallite_partionedLp0, & !< plastic velocity grad at start of homog inc - crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc ! crystallite_Li, & !< current intermediate velocitiy grad (end of converged time step) crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0, & !< intermediate velocity grad at start of homog inc - crystallite_subLi0, & !< intermediate velocity grad at start of crystallite inc ! crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc @@ -183,7 +181,6 @@ subroutine crystallite_init crystallite_Li,crystallite_Lp, & crystallite_subF,crystallite_subF0, & crystallite_subFp0,crystallite_subFi0, & - crystallite_subLi0,crystallite_subLp0, & source = crystallite_partionedF) allocate(crystallite_dPdF(3,3,3,3,cMax,iMax,eMax),source=0.0_pReal) @@ -326,34 +323,17 @@ function crystallite_stress() startIP, endIP, & s logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains - todo = .false. + real(pReal), dimension(:,:,:,:,:), allocatable :: & + subLp0,& !< plastic velocity grad at start of crystallite inc + subLi0 !< intermediate velocity grad at start of crystallite inc + + + todo = .false. + + subLp0 = crystallite_partionedLp0 + subLi0 = crystallite_partionedLi0 + -#ifdef DEBUG - if (debugCrystallite%selective & - .and. FEsolving_execElem(1) <= debugCrystallite%element & - .and. debugCrystallite%element <= FEsolving_execElem(2)) then - print'(/,a,i8,1x,i2,1x,i3)', '<< CRYST stress >> boundary and initial values at el ip ipc ', & - debugCrystallite%element,debugCrystallite%ip, debugCrystallite%grain - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F ', & - transpose(crystallite_partionedF(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> F0 ', & - transpose(crystallite_partionedF0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fp0', & - transpose(crystallite_partionedFp0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Fi0', & - transpose(crystallite_partionedFi0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Lp0', & - transpose(crystallite_partionedLp0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - print'(a,/,3(12x,3(f14.9,1x)/))', '<< CRYST stress >> Li0', & - transpose(crystallite_partionedLi0(1:3,1:3,debugCrystallite%grain, & - debugCrystallite%ip,debugCrystallite%element)) - endif -#endif !-------------------------------------------------------------------------------------------------- ! initialize to starting condition @@ -370,9 +350,7 @@ function crystallite_stress() sourceState(material_phaseAt(c,e))%p(s)%partionedState0(:,material_phaseMemberAt(c,i,e)) enddo crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_partionedFp0(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_partionedLp0(1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_partionedFi0(1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_partionedLi0(1:3,1:3,c,i,e) crystallite_subF0(1:3,1:3,c,i,e) = crystallite_partionedF0(1:3,1:3,c,i,e) crystallite_subFrac(c,i,e) = 0.0_pReal crystallite_subStep(c,i,e) = 1.0_pReal/num%subStepSizeCryst @@ -415,8 +393,8 @@ function crystallite_stress() todo(c,i,e) = crystallite_subStep(c,i,e) > 0.0_pReal ! still time left to integrate on? if (todo(c,i,e)) then crystallite_subF0 (1:3,1:3,c,i,e) = crystallite_subF(1:3,1:3,c,i,e) - crystallite_subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) - crystallite_subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) + subLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) + subLi0(1:3,1:3,c,i,e) = crystallite_Li (1:3,1:3,c,i,e) crystallite_subFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) crystallite_subFi0(1:3,1:3,c,i,e) = crystallite_Fi (1:3,1:3,c,i,e) plasticState( material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) & @@ -435,8 +413,8 @@ function crystallite_stress() crystallite_Fi (1:3,1:3,c,i,e) = crystallite_subFi0(1:3,1:3,c,i,e) crystallite_S (1:3,1:3,c,i,e) = crystallite_S0 (1:3,1:3,c,i,e) if (crystallite_subStep(c,i,e) < 1.0_pReal) then ! actual (not initial) cutback - crystallite_Lp (1:3,1:3,c,i,e) = crystallite_subLp0(1:3,1:3,c,i,e) - crystallite_Li (1:3,1:3,c,i,e) = crystallite_subLi0(1:3,1:3,c,i,e) + crystallite_Lp (1:3,1:3,c,i,e) = subLp0(1:3,1:3,c,i,e) + crystallite_Li (1:3,1:3,c,i,e) = subLi0(1:3,1:3,c,i,e) endif plasticState (material_phaseAt(c,e))%state( :,material_phaseMemberAt(c,i,e)) & = plasticState(material_phaseAt(c,e))%subState0(:,material_phaseMemberAt(c,i,e)) From 54e4943353829708cf09d97ac1e3d085a75d8c81 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 30 Sep 2020 07:10:15 +0200 Subject: [PATCH 29/39] get rid of shell scripts --- processing/pre/geom_fromMinimalSurface.py | 21 ++------------ python/damask/_geom.py | 35 +++++++++++++++++++++++ python/tests/test_Geom.py | 18 ++++++++++++ 3 files changed, 55 insertions(+), 19 deletions(-) diff --git a/processing/pre/geom_fromMinimalSurface.py b/processing/pre/geom_fromMinimalSurface.py index 83d1a3684..bf0e2b754 100755 --- a/processing/pre/geom_fromMinimalSurface.py +++ b/processing/pre/geom_fromMinimalSurface.py @@ -4,8 +4,6 @@ import os import sys from optparse import OptionParser -import numpy as np - import damask @@ -15,13 +13,6 @@ scriptID = ' '.join([scriptName,damask.version]) minimal_surfaces = ['primitive','gyroid','diamond'] -surface = { - 'primitive': lambda x,y,z: np.cos(x)+np.cos(y)+np.cos(z), - 'gyroid': lambda x,y,z: np.sin(x)*np.cos(y)+np.sin(y)*np.cos(z)+np.cos(x)*np.sin(z), - 'diamond': lambda x,y,z: np.cos(x-y)*np.cos(z)+np.sin(x+y)*np.sin(z), - } - - # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -71,16 +62,8 @@ parser.set_defaults(type = minimal_surfaces[0], name = None if filename == [] else filename[0] damask.util.report(scriptName,name) -x,y,z = np.meshgrid(options.periods*2.0*np.pi*(np.arange(options.grid[0])+0.5)/options.grid[0], - options.periods*2.0*np.pi*(np.arange(options.grid[1])+0.5)/options.grid[1], - options.periods*2.0*np.pi*(np.arange(options.grid[2])+0.5)/options.grid[2], - indexing='xy',sparse=True) - -microstructure = np.where(options.threshold < surface[options.type](x,y,z), - options.microstructure[1],options.microstructure[0]) - -geom=damask.Geom(microstructure,options.size, - comments=[scriptID + ' ' + ' '.join(sys.argv[1:])]) +geom=damask.Geom.from_minimal_surface(options.grid,options.size,options.type,options.threshold, + options.periods,options.microstructure) damask.util.croak(geom) geom.save_ASCII(sys.stdout if name is None else name,compress=False) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 6a7bc3c0e..aaa534df0 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -292,6 +292,41 @@ class Geom: ) + @staticmethod + def from_minimal_surface(grid,size,surface,threshold=0.0,periods=1,materials=(1,2)): + """ + Generate geometry from definition of minimal surface. + + Parameters + ---------- + grid : int numpy.ndarray of shape (3) + Number of grid points in x,y,z direction. + size : list or numpy.ndarray of shape (3) + Physical size of the geometry in meter. + surface : {'primitive', 'gyroid', 'diamond'} + Type of the minimal surface. + threshold : float, optional. + Threshold of the minimal surface. Defaults to 0.0. + periods : integer, optional. + Number of periods per unit cell. Defaults to 1. + materials : (int, int), optional + Material IDs. Defaults to (1,2). + + """ + s = {'primitive': lambda x,y,z: np.cos(x)+np.cos(y)+np.cos(z), + 'gyroid': lambda x,y,z: np.sin(x)*np.cos(y)+np.sin(y)*np.cos(z)+np.cos(x)*np.sin(z), + 'diamond': lambda x,y,z: np.cos(x-y)*np.cos(z)+np.sin(x+y)*np.sin(z), + } + x,y,z = np.meshgrid(periods*2.0*np.pi*(np.arange(grid[0])+0.5)/grid[0], + periods*2.0*np.pi*(np.arange(grid[1])+0.5)/grid[1], + periods*2.0*np.pi*(np.arange(grid[2])+0.5)/grid[2], + indexing='ij',sparse=True) + return Geom(material = np.where(threshold < s[surface](x,y,z),materials[1],materials[0]), + size = size, + comments = util.execution_stamp('Geom','from_minimal_surface'), + ) + + def save_ASCII(self,fname,compress=None): """ Writes a geom file. diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 21eb38cdd..27829ac34 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -360,3 +360,21 @@ class TestGeom: elif approach == 'Voronoi': geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5) assert np.all(geom.material == material) + + + @pytest.mark.parametrize('surface',['primitive','gyroid','diamond']) + def test_minimal_surface_basic_properties(self,surface): + grid = np.random.randint(60,100,3) + size = np.ones(3)+np.random.rand(3) + threshold = np.random.rand()-.5 + periods = np.random.randint(2)+1 + materials = np.random.randint(0,40,2) + geom = Geom.from_minimal_surface(grid,size,surface,threshold,periods,materials) + assert geom.material.max() == materials.max() and geom.material.min() == materials.min() \ + and (geom.size == size).all() and (geom.grid == grid).all() + + @pytest.mark.parametrize('surface',['primitive','gyroid','diamond']) + def test_minimal_surface_volume(self,surface): + grid = np.ones(3,dtype='i')*64 + geom = Geom.from_minimal_surface(grid,np.ones(3),surface) + assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5) From ff0fdfa0aa084334739d9c4e7ab768a1764aff9f Mon Sep 17 00:00:00 2001 From: Test User Date: Wed, 30 Sep 2020 09:00:48 +0200 Subject: [PATCH 30/39] [skip ci] updated version information after successful test of v3.0.0-alpha-389-ge1d459aab --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index f11394f7c..5e3f5a21a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-381-g63f2419e9 +v3.0.0-alpha-389-ge1d459aab From a59e64a8e459b779a439edb0fc4022084d1de1d5 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 30 Sep 2020 17:27:49 -0400 Subject: [PATCH 31/39] renamed TPMS and added more from additional references --- python/damask/_geom.py | 154 +++++++++++++++++++++++++------------- python/tests/test_Geom.py | 36 +++++++-- 2 files changed, 134 insertions(+), 56 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index aaa534df0..e760e40c3 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -120,6 +120,29 @@ class Geom: return np.unique(self.material).size + @staticmethod + def load(fname): + """ + Read a VTK rectilinear grid. + + Parameters + ---------- + fname : str or or pathlib.Path + Geometry file to read. + Valid extension is .vtr, it will be appended if not given. + + """ + v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr') + comments = v.get_comments() + grid = np.array(v.vtk_data.GetDimensions())-1 + bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T + + return Geom(material = v.get('material').reshape(grid,order='F'), + size = bbox[1] - bbox[0], + origin = bbox[0], + comments=comments) + + @staticmethod def load_ASCII(fname): """ @@ -184,29 +207,6 @@ class Geom: return Geom(material.reshape(grid,order='F'),size,origin,comments) - @staticmethod - def load(fname): - """ - Read a VTK rectilinear grid. - - Parameters - ---------- - fname : str or or pathlib.Path - Geometry file to read. - Valid extension is .vtr, it will be appended if not given. - - """ - v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr') - comments = v.get_comments() - grid = np.array(v.vtk_data.GetDimensions())-1 - bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T - - return Geom(material = v.get('material').reshape(grid,order='F'), - size = bbox[1] - bbox[0], - origin = bbox[0], - comments=comments) - - @staticmethod def _find_closest_seed(seeds, weights, point): return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) @@ -292,10 +292,52 @@ class Geom: ) + _minimal_surface = \ + {'Schwarz P': lambda x,y,z: np.cos(x) + np.cos(y) + np.cos(z), + 'Double Primitive': lambda x,y,z: ( 0.5 * (np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x)) + + 0.2 * (np.cos(2*x) + np.cos(2*y) + np.cos(2*z)) ), + 'Schwarz D': lambda x,y,z: ( np.sin(x)*np.sin(y)*np.sin(z) + + np.sin(x)*np.cos(y)*np.cos(z) + + np.cos(x)*np.cos(y)*np.sin(z) + + np.cos(x)*np.sin(y)*np.cos(z) ), + 'Complementary D': lambda x,y,z: ( np.cos(3*x+y)*np.cos(z) - np.sin(3*x-y)*np.sin(z) + np.cos(x+3*y)*np.cos(z) + + np.sin(x-3*y)*np.sin(z) + np.cos(x-y)*np.cos(3*z) - np.sin(x+y)*np.sin(3*z) ), + 'Double Diamond': lambda x,y,z: 0.5 * (np.sin(x)*np.sin(y) + + np.sin(y)*np.sin(z) + + np.sin(z)*np.sin(x) + + np.cos(x) * np.cos(y) * np.cos(z) ), + 'Dprime': lambda x,y,z: 0.5 * ( np.cos(x)*np.cos(y)*np.cos(z) + + np.cos(x)*np.sin(y)*np.sin(z) + + np.sin(x)*np.cos(y)*np.sin(z) + + np.sin(x)*np.sin(y)*np.cos(z) + - np.sin(2*x)*np.sin(2*y) + - np.sin(2*y)*np.sin(2*z) + - np.sin(2*z)*np.sin(2*x) ) - 0.2, + 'Gyroid': lambda x,y,z: np.cos(x)*np.sin(y) + np.cos(y)*np.sin(z) + np.cos(z)*np.sin(x), + 'Gprime': lambda x,y,z : ( np.sin(2*x)*np.cos(y)*np.sin(z) + + np.sin(2*y)*np.cos(z)*np.sin(x) + + np.sin(2*z)*np.cos(x)*np.sin(y) ) + 0.32, + 'Karcher K': lambda x,y,z: ( 0.3 * ( np.cos(x) + np.cos(y) + np.cos(z) + + np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x) ) + - 0.4 * ( np.cos(2*x) + np.cos(2*y) + np.cos(2*z) ) ) + 0.2, + 'Lidinoid': lambda x,y,z: 0.5 * ( np.sin(2*x)*np.cos(y)*np.sin(z) + + np.sin(2*y)*np.cos(z)*np.sin(x) + + np.sin(2*z)*np.cos(x)*np.sin(y) + - np.cos(2*x)*np.cos(2*y) + - np.cos(2*y)*np.cos(2*z) + - np.cos(2*z)*np.cos(2*x) ) + 0.15, + 'Neovius': lambda x,y,z: ( 3 * (np.cos(x)+np.cos(y)+np.cos(z)) + + 4 * np.cos(x)*np.cos(y)*np.cos(z) ), + 'Fisher-Koch S': lambda x,y,z: ( np.cos(2*x)*np.sin( y)*np.cos( z) + + np.cos( x)*np.cos(2*y)*np.sin( z) + + np.sin( x)*np.cos( y)*np.cos(2*z) ), + } + + @staticmethod def from_minimal_surface(grid,size,surface,threshold=0.0,periods=1,materials=(1,2)): """ - Generate geometry from definition of minimal surface. + Generate geometry from definition of triply periodic minimal surface. Parameters ---------- @@ -303,7 +345,7 @@ class Geom: Number of grid points in x,y,z direction. size : list or numpy.ndarray of shape (3) Physical size of the geometry in meter. - surface : {'primitive', 'gyroid', 'diamond'} + surface : {'primitive', 'gyroid', 'lidinoid', 'neovius', 'diamond', 'doublediamond'} Type of the minimal surface. threshold : float, optional. Threshold of the minimal surface. Defaults to 0.0. @@ -312,21 +354,53 @@ class Geom: materials : (int, int), optional Material IDs. Defaults to (1,2). + References + ---------- + Surface curvature in triply-periodic minimal surface architectures as + a distinct design parameter in preparing advanced tissue engineering scaffolds + Sébastien B G Blanquer, Maike Werner, Markus Hannula, Shahriar Sharifi, + Guillaume P R Lajoinie, David Eglin, Jari Hyttinen, André A Poot, and Dirk W Grijpma + 10.1088/1758-5090/aa6553 + + Triply Periodic Bicontinuous Cubic Microdomain Morphologies by Symmetries + Meinhard Wohlgemuth, Nataliya Yufa, James Hoffman, and Edwin L. Thomas + 10.1021/ma0019499 + + Minisurf – A minimal surface generator for finite element modeling and additive manufacturing + Meng-Ting Hsieh, Lorenzo Valdevit + 10.1016/j.simpa.2020.100026 + """ - s = {'primitive': lambda x,y,z: np.cos(x)+np.cos(y)+np.cos(z), - 'gyroid': lambda x,y,z: np.sin(x)*np.cos(y)+np.sin(y)*np.cos(z)+np.cos(x)*np.sin(z), - 'diamond': lambda x,y,z: np.cos(x-y)*np.cos(z)+np.sin(x+y)*np.sin(z), - } x,y,z = np.meshgrid(periods*2.0*np.pi*(np.arange(grid[0])+0.5)/grid[0], periods*2.0*np.pi*(np.arange(grid[1])+0.5)/grid[1], periods*2.0*np.pi*(np.arange(grid[2])+0.5)/grid[2], indexing='ij',sparse=True) - return Geom(material = np.where(threshold < s[surface](x,y,z),materials[1],materials[0]), + return Geom(material = np.where(threshold < Geom._minimal_surface[surface](x,y,z),materials[1],materials[0]), size = size, comments = util.execution_stamp('Geom','from_minimal_surface'), ) + def save(self,fname,compress=True): + """ + Generates vtk rectilinear grid. + + Parameters + ---------- + fname : str, optional + Filename to write. If no file is given, a string is returned. + Valid extension is .vtr, it will be appended if not given. + compress : bool, optional + Compress with zlib algorithm. Defaults to True. + + """ + v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) + v.add(self.material.flatten(order='F'),'material') + v.add_comments(self.comments) + + v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) + + def save_ASCII(self,fname,compress=None): """ Writes a geom file. @@ -399,26 +473,6 @@ class Geom: f.write(f'{reps} of {former}\n') - def save(self,fname,compress=True): - """ - Generates vtk rectilinear grid. - - Parameters - ---------- - fname : str, optional - Filename to write. If no file is given, a string is returned. - Valid extension is .vtr, it will be appended if not given. - compress : bool, optional - Compress with zlib algorithm. Defaults to True. - - """ - v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) - v.add(self.material.flatten(order='F'),'material') - v.add_comments(self.comments) - - v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) - - def show(self): """Show on screen.""" v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 27829ac34..884ae72e7 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -362,7 +362,19 @@ class TestGeom: assert np.all(geom.material == material) - @pytest.mark.parametrize('surface',['primitive','gyroid','diamond']) + @pytest.mark.parametrize('surface',['Schwarz P', + 'Double Primitive', + 'Schwarz D', + 'Complementary D', + 'Double Diamond', + 'Dprime', + 'Gyroid', + 'Gprime', + 'Karcher K', + 'Lidinoid', + 'Neovius', + 'Fisher-Koch S', + ]) def test_minimal_surface_basic_properties(self,surface): grid = np.random.randint(60,100,3) size = np.ones(3)+np.random.rand(3) @@ -373,8 +385,20 @@ class TestGeom: assert geom.material.max() == materials.max() and geom.material.min() == materials.min() \ and (geom.size == size).all() and (geom.grid == grid).all() - @pytest.mark.parametrize('surface',['primitive','gyroid','diamond']) - def test_minimal_surface_volume(self,surface): - grid = np.ones(3,dtype='i')*64 - geom = Geom.from_minimal_surface(grid,np.ones(3),surface) - assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5) + @pytest.mark.parametrize('surface,threshold',[('Schwarz P',0), + ('Double Primitive',-1./6.), + ('Schwarz D',0), + ('Complementary D',0), + ('Double Diamond',-0.133), + ('Dprime',-0.0395), + ('Gyroid',0), + ('Gprime',0.22913), + ('Karcher K',0.17045), + ('Lidinoid',0.14455), + ('Neovius',0), + ('Fisher-Koch S',0), + ]) + def test_minimal_surface_volume(self,surface,threshold): + grid = np.ones(3,dtype=int)*64 + geom = Geom.from_minimal_surface(grid,np.ones(3),surface,threshold) + assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5,rtol=1e-3) From f2048087954aaaf33f15f946e67653342e7d4c0d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 30 Sep 2020 17:45:02 -0400 Subject: [PATCH 32/39] adapted geom_fromMinimalSurface.py to new TPMS names --- processing/pre/geom_fromMinimalSurface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/pre/geom_fromMinimalSurface.py b/processing/pre/geom_fromMinimalSurface.py index bf0e2b754..2b40940b7 100755 --- a/processing/pre/geom_fromMinimalSurface.py +++ b/processing/pre/geom_fromMinimalSurface.py @@ -11,7 +11,7 @@ scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -minimal_surfaces = ['primitive','gyroid','diamond'] +minimal_surfaces = list(damask.Geom._minimal_surface.keys()) # -------------------------------------------------------------------- # MAIN From 65fe3f821ac1f02266009d6579c929617dd6d7fa Mon Sep 17 00:00:00 2001 From: Test User Date: Thu, 1 Oct 2020 01:42:33 +0200 Subject: [PATCH 33/39] [skip ci] updated version information after successful test of v3.0.0-alpha-400-ga3674d931 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 5e3f5a21a..9d7d6ae9a 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-389-ge1d459aab +v3.0.0-alpha-400-ga3674d931 From b29f22f51398b7c36608ab5b693239e00c19cb41 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 1 Oct 2020 09:25:32 +0200 Subject: [PATCH 34/39] documenting the actually available TMPS --- python/damask/_geom.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index e760e40c3..2ccbb1988 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -345,14 +345,30 @@ class Geom: Number of grid points in x,y,z direction. size : list or numpy.ndarray of shape (3) Physical size of the geometry in meter. - surface : {'primitive', 'gyroid', 'lidinoid', 'neovius', 'diamond', 'doublediamond'} - Type of the minimal surface. + surface : str + Type of the minimal surface. See notes for details. threshold : float, optional. Threshold of the minimal surface. Defaults to 0.0. periods : integer, optional. Number of periods per unit cell. Defaults to 1. materials : (int, int), optional Material IDs. Defaults to (1,2). + + Notes + ----- + The following triply-periodic minimal surfaces are implemented: + - Schwarz P + - Double Primitive + - Schwarz D + - Complementary D + - Double Diamond + - Dprime + - Gyroid + - Gprime + - Karcher K + - Lidinoid + - Neovius + - Fisher-Koch S References ---------- From 776901cb81f569b801f702f1abaf5d24d92e37fc Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 1 Oct 2020 12:43:05 +0200 Subject: [PATCH 35/39] microstructure --> material --- PRIVATE | 2 +- .../SpectralMethod/Polycrystal/material.yaml | 207 +++++++++--------- src/IO.f90 | 10 +- src/discretization.f90 | 10 +- src/material.f90 | 32 +-- src/mesh/discretization_mesh.f90 | 8 +- 6 files changed, 135 insertions(+), 134 deletions(-) diff --git a/PRIVATE b/PRIVATE index 25ce39dd0..8f27fc91c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 25ce39dd0f5bf49dc5c2bec20767d93d2b76d353 +Subproject commit 8f27fc91ca757a1dfdfd04892708af7e94941ef9 diff --git a/examples/SpectralMethod/Polycrystal/material.yaml b/examples/SpectralMethod/Polycrystal/material.yaml index 16c6042a6..99ff0440c 100644 --- a/examples/SpectralMethod/Polycrystal/material.yaml +++ b/examples/SpectralMethod/Polycrystal/material.yaml @@ -1,107 +1,109 @@ homogenization: SX: mech: {type: none} -microstructure: -- constituents: - - fraction: 1.0 - orientation: [1.0, 0.0, 0.0, 0.0] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543] - phase: Aluminum - homogenization: SX -- constituents: - - fraction: 1.0 - orientation: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341] - phase: Aluminum - homogenization: SX + +material: + - constituents: + - fraction: 1.0 + orientation: [1.0, 0.0, 0.0, 0.0] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543] + phase: Aluminum + homogenization: SX + - constituents: + - fraction: 1.0 + orientation: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341] + phase: Aluminum + homogenization: SX + phase: Aluminum: elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke} @@ -117,7 +119,6 @@ phase: h_sl_sl: [1, 1, 1.4, 1.4, 1.4, 1.4] n_sl: 20 output: [xi_sl] + type: phenopowerlaw xi_0_sl: [31e6] xi_inf_sl: [63e6] - type: phenopowerlaw - diff --git a/src/IO.f90 b/src/IO.f90 index 0542e7a62..4a442cb06 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -427,20 +427,20 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (146) msg = 'number of values does not match' case (148) - msg = 'Nconstituents mismatch between homogenization and microstructure' + msg = 'Nconstituents mismatch between homogenization and material' !-------------------------------------------------------------------------------------------------- ! material error messages and related messages in mesh case (150) msg = 'index out of bounds' case (151) - msg = 'microstructure has no constituents' + msg = 'material has no constituents' case (153) msg = 'sum of phase fractions differs from 1' case (155) - msg = 'microstructure index out of bounds' + msg = 'material index out of bounds' case (180) - msg = 'missing/invalid microstructure definition via State Variable 2' + msg = 'missing/invalid material definition via State Variable 2' case (190) msg = 'unknown element type:' case (191) @@ -526,7 +526,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (842) msg = 'incomplete information in grid mesh header' case (843) - msg = 'microstructure count mismatch' + msg = 'material count mismatch' case (844) msg = 'invalid VTR file' case (846) diff --git a/src/discretization.f90 b/src/discretization.f90 index 3bb905cf3..e6e53fcf4 100644 --- a/src/discretization.f90 +++ b/src/discretization.f90 @@ -15,7 +15,7 @@ module discretization discretization_nElem integer, public, protected, dimension(:), allocatable :: & - discretization_microstructureAt + discretization_materialAt real(pReal), public, protected, dimension(:,:), allocatable :: & discretization_IPcoords0, & @@ -37,12 +37,12 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief stores the relevant information in globally accesible variables !-------------------------------------------------------------------------------------------------- -subroutine discretization_init(microstructureAt,& +subroutine discretization_init(materialAt,& IPcoords0,NodeCoords0,& sharedNodesBegin) integer, dimension(:), intent(in) :: & - microstructureAt + materialAt real(pReal), dimension(:,:), intent(in) :: & IPcoords0, & NodeCoords0 @@ -51,10 +51,10 @@ subroutine discretization_init(microstructureAt,& print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6) - discretization_nElem = size(microstructureAt,1) + discretization_nElem = size(materialAt,1) discretization_nIP = size(IPcoords0,2)/discretization_nElem - discretization_microstructureAt = microstructureAt + discretization_materialAt = materialAt discretization_IPcoords0 = IPcoords0 discretization_IPcoords = IPcoords0 diff --git a/src/material.f90 b/src/material.f90 index 6688a0ae5..3a2483246 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -97,7 +97,7 @@ module material material_orientation0 !< initial orientation of each grain,IP,element integer, dimension(:), allocatable, private :: & - microstructure_Nconstituents !< number of constituents in each microstructure + material_Nconstituents !< number of constituents in each material @@ -317,12 +317,12 @@ end subroutine material_parseHomogenization !-------------------------------------------------------------------------------------------------- -!> @brief parses the microstructure part in the material configuration file +!> @brief parses the material part in the material configuration file !-------------------------------------------------------------------------------------------------- subroutine material_parseMicrostructure - class(tNode), pointer :: microstructures, & !> list of microstructures - microstructure, & !> microstructure definition + class(tNode), pointer :: materials, & !> list of materials + material, & !> material definition constituents, & !> list of constituents constituent, & !> constituent definition phases, & @@ -341,17 +341,17 @@ subroutine material_parseMicrostructure c, & maxNconstituents - microstructures => config_material%get('microstructure') - if(any(discretization_microstructureAt > microstructures%length)) & - call IO_error(155,ext_msg='More microstructures requested than found in material.yaml') + materials => config_material%get('material') + if(any(discretization_materialAt > materials%length)) & + call IO_error(155,ext_msg='More materials requested than found in material.yaml') - allocate(microstructure_Nconstituents(microstructures%length),source=0) - do m = 1, microstructures%length - microstructure => microstructures%get(m) - constituents => microstructure%get('constituents') - microstructure_Nconstituents(m) = constituents%length + allocate(material_Nconstituents(materials%length),source=0) + do m = 1, materials%length + material => materials%get(m) + constituents => material%get('constituents') + material_Nconstituents(m) = constituents%length enddo - maxNconstituents = maxval(microstructure_Nconstituents) + maxNconstituents = maxval(material_Nconstituents) allocate(material_homogenizationAt(discretization_nElem),source=0) allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0) @@ -366,10 +366,10 @@ subroutine material_parseMicrostructure allocate(counterHomogenization(homogenization%length),source=0) do e = 1, discretization_nElem - microstructure => microstructures%get(discretization_microstructureAt(e)) - constituents => microstructure%get('constituents') + material => materials%get(discretization_materialAt(e)) + constituents => material%get('constituents') - material_homogenizationAt(e) = homogenization%getIndex(microstructure%get_asString('homogenization')) + material_homogenizationAt(e) = homogenization%getIndex(material%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)) diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index b57a4b332..bc96951a1 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -78,7 +78,7 @@ subroutine discretization_mesh_init(restart) IS :: faceSetIS PetscErrorCode :: ierr integer, dimension(:), allocatable :: & - microstructureAt + materialAt class(tNode), pointer :: & num_mesh integer :: integrationOrder !< order of quadrature rule required @@ -164,9 +164,9 @@ subroutine discretization_mesh_init(restart) call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p) call mesh_FEM_build_ipVolumes(dimPlex) - allocate(microstructureAt(mesh_NcpElems)) + allocate(materialAt(mesh_NcpElems)) do j = 1, mesh_NcpElems - call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr) + call DMGetLabelValue(geomMesh,'material',j-1,materialAt(j),ierr) CHKERRQ(ierr) end do @@ -178,7 +178,7 @@ subroutine discretization_mesh_init(restart) allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal) - call discretization_init(microstructureAt,& + call discretization_init(materialAt,& reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), & mesh_node0) From 33b0181286a5b51b48f026ee6356c7444c02a097 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 1 Oct 2020 14:01:50 +0200 Subject: [PATCH 36/39] orientation --> O --- PRIVATE | 2 +- .../SpectralMethod/Polycrystal/material.yaml | 40 +++++++++---------- src/material.f90 | 4 +- 3 files changed, 23 insertions(+), 23 deletions(-) diff --git a/PRIVATE b/PRIVATE index 8f27fc91c..28e7fb09f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 8f27fc91ca757a1dfdfd04892708af7e94941ef9 +Subproject commit 28e7fb09f06ec324841eac1670ebdcbd8aab4ef3 diff --git a/examples/SpectralMethod/Polycrystal/material.yaml b/examples/SpectralMethod/Polycrystal/material.yaml index 99ff0440c..394d16ed8 100644 --- a/examples/SpectralMethod/Polycrystal/material.yaml +++ b/examples/SpectralMethod/Polycrystal/material.yaml @@ -5,102 +5,102 @@ homogenization: material: - constituents: - fraction: 1.0 - orientation: [1.0, 0.0, 0.0, 0.0] + O: [1.0, 0.0, 0.0, 0.0] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434] + O: [0.7936696712125002, -0.28765777461664166, -0.3436487135089419, 0.4113964260949434] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945] + O: [0.3986143167493579, -0.7014883552495493, 0.2154871765709027, 0.5500781677772945] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106] + O: [0.28645844315788244, -0.022571491243423537, -0.467933059311115, -0.8357456192708106] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676] + O: [0.33012772942625784, -0.6781865350268957, 0.6494525351030648, 0.09638521992649676] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265] + O: [0.43596817439583935, -0.5982537129781701, 0.046599032277502436, 0.6707106499919265] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854] + O: [0.169734823419553, -0.699615227367322, -0.6059581215838098, -0.33844257746495854] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175] + O: [0.9698864809294915, 0.1729052643205874, -0.15948307917616958, 0.06315956884687175] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738] + O: [0.46205660912967883, 0.3105054068891252, -0.617849551030653, 0.555294529545738] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815] + O: [0.4512443497461787, -0.7636045534540555, -0.04739348426715133, -0.45939142396805815] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289] + O: [0.2161856212656443, -0.6581450184826598, -0.5498086209601588, 0.4667112513346289] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541] + O: [0.8753220715350803, -0.4561599367657419, -0.13298279533852678, -0.08969369719975541] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101] + O: [0.11908260752431069, 0.18266024809834172, -0.7144822594012615, -0.664807992845101] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861] + O: [0.751104669484278, 0.5585633382623958, -0.34579336397009175, 0.06538900566860861] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363] + O: [0.08740438971703973, 0.8991264096610437, -0.4156704205935976, 0.10559485570696363] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859] + O: [0.5584325870096193, 0.6016408353068798, -0.14280340445801173, 0.5529814994483859] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182] + O: [0.4052725440888093, 0.25253073423599154, 0.5693263597910454, -0.669215876471182] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105] + O: [0.7570164606888676, 0.15265448024694664, -0.5998021466848317, 0.20942796551297105] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543] + O: [0.6987659297138081, -0.132172211261028, -0.19693254724422338, 0.6748883269678543] phase: Aluminum homogenization: SX - constituents: - fraction: 1.0 - orientation: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341] + O: [0.7729330445886478, 0.21682179052722322, -0.5207379472917645, 0.2905078484066341] phase: Aluminum homogenization: SX diff --git a/src/material.f90 b/src/material.f90 index 3a2483246..eece6626e 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -97,7 +97,7 @@ module material material_orientation0 !< initial orientation of each grain,IP,element integer, dimension(:), allocatable, private :: & - material_Nconstituents !< number of constituents in each material + material_Nconstituents !< number of constituents in each material @@ -385,7 +385,7 @@ subroutine material_parseMicrostructure counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) - call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('orientation',requiredSize=4)) + call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) enddo enddo From 26acdaf9ecba588316a43b06a703cdfd75bfa7cf Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 1 Oct 2020 16:07:50 +0200 Subject: [PATCH 37/39] missed out --- src/material.f90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/material.f90 b/src/material.f90 index eece6626e..81f46d64e 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -180,8 +180,8 @@ subroutine material_init(restart) material_name_homogenization(myHomog) = trim(adjustl(sectionName))//material_homogenization%getKey(myHomog) enddo - call material_parseMicrostructure - print*, 'Microstructure parsed' + call material_parseMaterial + print*, 'Material parsed' call material_parseHomogenization print*, 'Homogenization parsed' @@ -319,7 +319,7 @@ end subroutine material_parseHomogenization !-------------------------------------------------------------------------------------------------- !> @brief parses the material part in the material configuration file !-------------------------------------------------------------------------------------------------- -subroutine material_parseMicrostructure +subroutine material_parseMaterial class(tNode), pointer :: materials, & !> list of materials material, & !> material definition @@ -393,7 +393,7 @@ subroutine material_parseMicrostructure enddo -end subroutine material_parseMicrostructure +end subroutine material_parseMaterial end module material From 4a7701d86a89d630f18e6eecf9a4f8123600852c Mon Sep 17 00:00:00 2001 From: Test User Date: Fri, 2 Oct 2020 03:46:25 +0200 Subject: [PATCH 38/39] [skip ci] updated version information after successful test of v3.0.0-alpha-406-gca6f7f7a3 --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/VERSION b/VERSION index 9d7d6ae9a..8ec654de6 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-400-ga3674d931 +v3.0.0-alpha-406-gca6f7f7a3 From 58229b8851fac4195250e57ce01e06a78f1fe909 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 2 Oct 2020 13:26:11 +0200 Subject: [PATCH 39/39] relaxed test conditions even threshold between -.5 and +.5 can result in a single material. --- python/tests/test_Geom.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 884ae72e7..c2e476f8f 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -378,11 +378,11 @@ class TestGeom: def test_minimal_surface_basic_properties(self,surface): grid = np.random.randint(60,100,3) size = np.ones(3)+np.random.rand(3) - threshold = np.random.rand()-.5 + threshold = 2*np.random.rand()-1. periods = np.random.randint(2)+1 materials = np.random.randint(0,40,2) geom = Geom.from_minimal_surface(grid,size,surface,threshold,periods,materials) - assert geom.material.max() == materials.max() and geom.material.min() == materials.min() \ + assert set(geom.material.flatten()) | set(materials) == set(materials) \ and (geom.size == size).all() and (geom.grid == grid).all() @pytest.mark.parametrize('surface,threshold',[('Schwarz P',0),