From 872f6a9d90df1a4ffd49dea319456bff35a017af Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 13 May 2013 09:44:23 +0000 Subject: [PATCH] introduced possibility so specify deformation gradient aim at end of load case, rate will be calculated using difference between start of load case and aim. Needed for cyclic loading. Use keyword "f" for this behavior, don't use it as short name for freq any more! --- code/DAMASK_spectral_driver.f90 | 34 ++++++++++++----------- code/DAMASK_spectral_solverAL.f90 | 20 +++++++------ code/DAMASK_spectral_solverBasic.f90 | 9 ++++-- code/DAMASK_spectral_solverBasicPETSc.f90 | 14 +++++++--- 4 files changed, 46 insertions(+), 31 deletions(-) diff --git a/code/DAMASK_spectral_driver.f90 b/code/DAMASK_spectral_driver.f90 index c325431d1..30b38d22c 100644 --- a/code/DAMASK_spectral_driver.f90 +++ b/code/DAMASK_spectral_driver.f90 @@ -108,10 +108,9 @@ program DAMASK_spectral_Driver integer(pInt), dimension(1_pInt + maxNchunks*2_pInt) :: positions ! this is longer than needed for geometry parsing integer(pInt) :: & - N_l = 0_pInt, & !< # of L specifiers found in load case file N_t = 0_pInt, & !< # of time indicators found in load case file N_n = 0_pInt, & !< # of increment specifiers found in load case file - N_Fdot = 0_pInt !< # of rate of F specifiers found in load case file + N_def = 0_pInt !< # of rate of deformation specifiers found in load case file character(len=1024) :: & line @@ -126,7 +125,8 @@ program DAMASK_spectral_Driver time = 0.0_pReal, & !< elapsed time time0 = 0.0_pReal, & !< begin of interval timeinc = 1.0_pReal, & !< current time interval - timeinc_old = 0.0_pReal ! !< previous time interval + timeIncOld = 0.0_pReal, & !< previous time interval + remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case logical :: & guess !< guess along former trajectory integer(pInt) :: & @@ -166,10 +166,8 @@ program DAMASK_spectral_Driver positions = IO_stringPos(line,maxNchunks) do i = 1_pInt, positions(1) ! reading compulsory parameters for loadcase select case (IO_lc(IO_stringValue(line,positions,i))) - case('l','velocitygrad','velgrad','velocitygradient') - N_l = N_l + 1_pInt - case('fdot','dotf') - N_Fdot = N_Fdot + 1_pInt + case('l','velocitygrad','velgrad','velocitygradient','fdot','dotf','f') + N_def = N_def + 1_pInt case('t','time','delta') N_t = N_t + 1_pInt case('n','incs','increments','steps','logincs','logincrements','logsteps') @@ -178,7 +176,7 @@ program DAMASK_spectral_Driver enddo ! count all identifiers to allocate memory and do sanity check enddo -100 if ((N_l + N_Fdot /= N_n) .or. (N_n /= N_t)) & ! sanity check +100 if ((N_def /= N_n) .or. (N_n /= N_t)) & ! sanity check call IO_error(error_ID=837_pInt,ext_msg = trim(loadCaseFile)) ! error message for incomplete loadcase allocate (loadCases(N_n)) ! array of load cases loadCases%P%myType='p' @@ -193,11 +191,13 @@ program DAMASK_spectral_Driver positions = IO_stringPos(line,maxNchunks) do i = 1_pInt, positions(1) select case (IO_lc(IO_stringValue(line,positions,i))) - case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient') ! assign values for the deformation BC matrix + case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient','f') ! assign values for the deformation BC matrix temp_valueVector = 0.0_pReal if (IO_lc(IO_stringValue(line,positions,i)) == 'fdot'.or. & ! in case of Fdot, set type to fdot IO_lc(IO_stringValue(line,positions,i)) == 'dotf') then loadCases(currentLoadCase)%deformation%myType = 'fdot' + else if (IO_lc(IO_stringValue(line,positions,i)) == 'f') then + loadCases(currentLoadCase)%deformation%myType = 'f' else loadCases(currentLoadCase)%deformation%myType = 'l' endif @@ -215,7 +215,7 @@ program DAMASK_spectral_Driver case('p','pk1','piolakirchhoff','stress', 's') temp_valueVector = 0.0_pReal do j = 1_pInt, 9_pInt - temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' ! true if not a * + temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' ! true if not an asterisk enddo do j = 1_pInt,9_pInt if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j) ! read value where applicable @@ -282,6 +282,8 @@ program DAMASK_spectral_Driver errorID = 832_pInt ! each row should be either fully or not at all defined enddo write(6,'(2x,a)') 'velocity gradient:' + else if (loadCases(currentLoadCase)%deformation%myType=='f') then + write(6,'(2x,a)') 'deformation gradient at end of load case:' else write(6,'(2x,a)') 'deformation gradient rate:' endif @@ -390,7 +392,7 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! forwarding time - timeinc_old = timeinc + timeIncOld = timeinc if (loadCases(currentLoadCase)%logscale == 0_pInt) then ! linear scale timeinc = loadCases(currentLoadCase)%time/loadCases(currentLoadCase)%incs ! only valid for given linear time scale. will be overwritten later in case loglinear scale is used else @@ -418,7 +420,7 @@ program DAMASK_spectral_Driver subIncLooping: do while (stepFraction/subStepFactor**cutBackLevel <1_pInt) time = time + timeinc ! forward time stepFraction = stepFraction + 1_pInt - + remainingLoadCaseTime = time0 - time + loadCases(currentLoadCase)%time + timeInc !-------------------------------------------------------------------------------------------------- ! report begin of new increment write(6,'(/,a)') ' ###########################################################################' @@ -441,7 +443,7 @@ program DAMASK_spectral_Driver ! calculate solution case (DAMASK_spectral_SolverBasic_label) solres = basic_solution (& - incInfo, guess,timeinc,timeinc_old, & + incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & P_BC = loadCases(currentLoadCase)%P, & F_BC = loadCases(currentLoadCase)%deformation, & temperature_bc = loadCases(currentLoadCase)%temperature, & @@ -449,7 +451,7 @@ program DAMASK_spectral_Driver #ifdef PETSc case (DAMASK_spectral_SolverBasicPETSC_label) solres = BasicPETSC_solution (& - incInfo, guess,timeinc,timeinc_old, & + incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & P_BC = loadCases(currentLoadCase)%P, & F_BC = loadCases(currentLoadCase)%deformation, & temperature_bc = loadCases(currentLoadCase)%temperature, & @@ -457,7 +459,7 @@ program DAMASK_spectral_Driver case (DAMASK_spectral_SolverAL_label) solres = AL_solution (& - incInfo, guess,timeinc,timeinc_old, & + incInfo,guess,timeinc,timeIncOld,remainingLoadCaseTime, & P_BC = loadCases(currentLoadCase)%P, & F_BC = loadCases(currentLoadCase)%deformation, & temperature_bc = loadCases(currentLoadCase)%temperature, & @@ -475,7 +477,7 @@ program DAMASK_spectral_Driver stepFraction = (stepFraction - 1_pInt) * subStepFactor ! adjust to new denominator cutBackLevel = cutBackLevel + 1_pInt time = time - timeinc ! rewind time - timeinc_old = timeinc + timeIncOld = timeinc timeinc = timeinc/2.0_pReal elseif (solres%termIll) then ! material point model cannot find a solution if(regridMode > 0_pInt) call quit(-1_pInt*(lastRestartWritten+1_pInt)) ! regrid requested (mode 1 or 2) diff --git a/code/DAMASK_spectral_solverAL.f90 b/code/DAMASK_spectral_solverAL.f90 index 9b27de989..1e21caf62 100644 --- a/code/DAMASK_spectral_solverAL.f90 +++ b/code/DAMASK_spectral_solverAL.f90 @@ -43,8 +43,8 @@ module DAMASK_spectral_solverAL DAMASK_spectral_SolverAL_label = 'al' !-------------------------------------------------------------------------------------------------- -! derived types ToDo: use here the type definition for a full loadcase including mask - type tSolutionParams +! derived types + type tSolutionParams !< @ToDo: use here the type definition for a full loadcase including mask real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal) :: timeinc real(pReal) :: temperature @@ -261,7 +261,7 @@ end subroutine AL_init !> @brief solution for the AL scheme with internal iterations !-------------------------------------------------------------------------------------------------- type(tSolutionState) function & - AL_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) + AL_solution(incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC) use numerics, only: & update_gamma, & itmax @@ -290,11 +290,13 @@ use mesh, only: & implicit none #include #include + !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: & - timeinc, & - timeinc_old, & + real(pReal), intent(in) :: & + timeinc, & !< increment in time for current solution + timeinc_old, & !< increment in time of last increment + loadCaseTime, & !< remaining time of current load case temperature_bc logical, intent(in) :: & guess @@ -346,7 +348,7 @@ use mesh, only: & endif AL_solution%converged =.false. - if ( cutBack) then + if (cutBack) then F_aim = F_aim_lastInc F_tau= reshape(F_tau_lastInc,[9,grid(1),grid(2),grid(3)]) F = reshape(F_lastInc, [9,grid(1),grid(2),grid(3)]) @@ -357,8 +359,10 @@ use mesh, only: & ! calculate rate for aim if (F_BC%myType=='l') then ! calculate f_aimDot from given L and current F f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed + elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed f_aimDot = F_BC%maskFloat * F_BC%values + elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime endif if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim diff --git a/code/DAMASK_spectral_solverBasic.f90 b/code/DAMASK_spectral_solverBasic.f90 index efa98fb0e..328779890 100644 --- a/code/DAMASK_spectral_solverBasic.f90 +++ b/code/DAMASK_spectral_solverBasic.f90 @@ -164,8 +164,8 @@ end subroutine basic_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the basic scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function & - basic_solution(incInfo,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) +type(tSolutionState) function basic_solution(& + incInfo,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC) use numerics, only: & itmax, & itmin, & @@ -212,7 +212,8 @@ type(tSolutionState) function & ! input data for solution real(pReal), intent(in) :: & timeinc, & !< increment in time for current solution - timeinc_old !< increment in time of last increment + timeinc_old, & !< increment in time of last increment + loadCaseTime !< remaining time of current load case logical, intent(in) :: & guess !< if .false., assume homogeneous addon type(tBoundaryCondition), intent(in) :: & @@ -274,6 +275,8 @@ type(tSolutionState) function & f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed f_aimDot = F_BC%maskFloat * F_BC%values + elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = F_BC%maskFloat * (F_BC%values - F_aim)/loadCaseTime endif if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim diff --git a/code/DAMASK_spectral_solverBasicPETSc.f90 b/code/DAMASK_spectral_solverBasicPETSc.f90 index 70d9000f8..998abf3b6 100644 --- a/code/DAMASK_spectral_solverBasicPETSc.f90 +++ b/code/DAMASK_spectral_solverBasicPETSc.f90 @@ -231,8 +231,8 @@ end subroutine basicPETSc_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the Basic PETSC scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(tSolutionState) function & - basicPETSc_solution(incInfoIn,guess,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) +type(tSolutionState) function basicPETSc_solution( & + incInfoIn,guess,timeinc,timeinc_old,loadCaseTime,P_BC,F_BC,temperature_bc,rotation_BC) use numerics, only: & update_gamma, & itmax @@ -262,7 +262,11 @@ type(tSolutionState) function & #include !-------------------------------------------------------------------------------------------------- ! input data for solution - real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc + real(pReal), intent(in) :: & + timeinc, & !< increment in time for current solution + timeinc_old, & !< increment in time of last increment + loadCaseTime, & !< remaining time of current load case + temperature_bc logical, intent(in):: guess type(tBoundaryCondition), intent(in) :: P_BC,F_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC @@ -314,7 +318,9 @@ type(tSolutionState) function & f_aimDot = F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) elseif(F_BC%myType=='fdot') then ! f_aimDot is prescribed f_aimDot = F_BC%maskFloat * F_BC%values - endif + elseif(F_BC%myType=='f') then ! aim at end of load case is prescribed + f_aimDot = F_BC%maskFloat * (F_BC%values -F_aim)/loadCaseTime + endif if (guess) f_aimDot = f_aimDot + P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old F_aim_lastInc = F_aim