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,:,:,:)