diff --git a/code/DAMASK_spectral_Driver.f90 b/code/DAMASK_spectral_Driver.f90 index 77470184a..745a752de 100644 --- a/code/DAMASK_spectral_Driver.f90 +++ b/code/DAMASK_spectral_Driver.f90 @@ -49,6 +49,7 @@ program DAMASK_spectral_Driver restartInc use numerics, only: & + maxCutBack, & rotation_tol, & mySpectralSolver @@ -57,9 +58,10 @@ program DAMASK_spectral_Driver materialpoint_results use DAMASK_spectral_Utilities, only: & - boundaryCondition, & - solutionState, & - debugGeneral + tBoundaryCondition, & + tSolutionState, & + debugGeneral, & + cutBack use DAMASK_spectral_SolverBasic #ifdef PETSc @@ -68,10 +70,9 @@ program DAMASK_spectral_Driver #endif implicit none - - type loadCase + type tLoadCase real(pReal), dimension (3,3) :: rotation = math_I3 ! rotation of BC - type(boundaryCondition) :: P, & ! stress BC + type(tBoundaryCondition) :: P, & ! stress BC deformation ! deformation BC (Fdot or L) real(pReal) :: time = 0.0_pReal, & ! length of increment temperature = 300.0_pReal ! isothermal starting conditions @@ -80,7 +81,7 @@ program DAMASK_spectral_Driver restartfrequency = 0_pInt, & ! frequency of restart writes logscale = 0_pInt ! linear/logaritmic time inc flag logical :: followFormerTrajectory = .true. ! follow trajectory of former loadcase - end type loadCase + end type tLoadCase !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file @@ -88,17 +89,15 @@ program DAMASK_spectral_Driver logical, dimension(9) :: temp_maskVector !> temporarily from loadcase file when reading in tensors integer(pInt), parameter :: maxNchunksLoadcase = (1_pInt + 9_pInt)*3_pInt +& ! deformation, rotation, and stress (1_pInt + 1_pInt)*5_pInt +& ! time, (log)incs, temp, restartfrequency, and outputfrequency - 1_pInt, & ! dropguessing - maxNchunksGeom = 7_pInt, & ! 4 identifiers, 3 values - myUnit = 234_pInt + 1_pInt ! dropguessing integer(pInt), dimension(1_pInt + maxNchunksLoadcase*2_pInt) :: positions ! this is longer than needed for geometry parsing integer(pInt) :: & N_l = 0_pInt, & N_t = 0_pInt, & N_n = 0_pInt, & - N_Fdot = 0_pInt ! number of Fourier points - + N_Fdot = 0_pInt, & ! number of Fourier points + myUnit = 234_pInt character(len=1024) :: & line @@ -107,16 +106,15 @@ program DAMASK_spectral_Driver ! loop variables, convergence etc. real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal real(pReal) :: time = 0.0_pReal, time0 = 0.0_pReal, timeinc = 1.0_pReal, timeinc_old = 0.0_pReal ! elapsed time, begin of interval, time interval, previous time interval - real(pReal) :: guessmode - real(pReal), dimension(3,3) :: temp33_Real - integer(pInt) :: i, j, k, l, errorID + real(pReal) :: guessmode + integer(pInt) :: i, j, k, l, errorID, cutBackLevel = 0_pInt, stepFraction = 0_pInt integer(pInt) :: currentLoadcase = 0_pInt, inc, & totalIncsCounter = 0_pInt,& notConvergedCounter = 0_pInt, convergedCounter = 0_pInt character(len=6) :: loadcase_string - - type(loadCase), allocatable, dimension(:) :: loadCases - type(solutionState) solres + character(len=1024) :: incInfo + type(tLoadCase), allocatable, dimension(:) :: loadCases + type(tSolutionState) solres !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) @@ -126,7 +124,6 @@ program DAMASK_spectral_Driver write(6,'(a)') ' <<<+- DAMASK_spectral_Driver init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" - write(6,'(a)') '' !-------------------------------------------------------------------------------------------------- ! reading basic information from load case file and allocate data structure containing load cases call IO_open_file(myUnit,trim(loadCaseFile)) @@ -234,8 +231,7 @@ program DAMASK_spectral_Driver errorID = 0_pInt checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases) write (loadcase_string, '(i6)' ) currentLoadCase - - write(6,'(2x,a,i6)') 'load case: ', currentLoadCase + write(6,'(1x,a,i6)') 'load case: ', currentLoadCase if (.not. loadCases(currentLoadCase)%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory' if (loadCases(currentLoadCase)%deformation%myType=='l') then @@ -247,20 +243,20 @@ program DAMASK_spectral_Driver else write(6,'(2x,a)') 'deformation gradient rate:' endif - write (6,'(3(3(3x,f12.7,1x)/))',advance='no') merge(math_transpose33(loadCases(currentLoadCase)%deformation%values),& - reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%deformation%maskLogical)) - write (6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'stress / GPa:',& + write(6,'(3(3(3x,f12.7,1x)/))',advance='no') merge(math_transpose33(loadCases(currentLoadCase)%deformation%values),& + reshape(spread(huge(1.0_pReal),1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%deformation%maskLogical)) + write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'stress / GPa:',& 1e-9_pReal*merge(math_transpose33(loadCases(currentLoadCase)%P%values),& - reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%P%maskLogical)) + reshape(spread(huge(1.0_pReal),1,9),[ 3,3]),transpose(loadCases(currentLoadCase)%P%maskLogical)) if (any(loadCases(currentLoadCase)%rotation /= math_I3)) & - write (6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& + write(6,'(2x,a,/,3(3(3x,f12.7,1x)/))',advance='no') 'rotation of loadframe:',& math_transpose33(loadCases(currentLoadCase)%rotation) write(6,'(2x,a,f12.6)') 'temperature:', loadCases(currentLoadCase)%temperature write(6,'(2x,a,f12.6)') 'time: ', loadCases(currentLoadCase)%time write(6,'(2x,a,i5)') 'increments: ', loadCases(currentLoadCase)%incs write(6,'(2x,a,i5)') 'output frequency: ', loadCases(currentLoadCase)%outputfrequency - write(6,'(2x,a,i5)') 'restart frequency: ', loadCases(currentLoadCase)%restartfrequency - + write(6,'(2x,a,i5,/)') 'restart frequency: ', loadCases(currentLoadCase)%restartfrequency + if (any(loadCases(currentLoadCase)%P%maskLogical .eqv. loadCases(currentLoadCase)%deformation%maskLogical)) errorID = 831_pInt ! exclusive or masking only if (any(loadCases(currentLoadCase)%P%maskLogical .and. transpose(loadCases(currentLoadCase)%P%maskLogical) .and. & reshape([ .false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false.],[ 3,3]))) & @@ -281,8 +277,8 @@ program DAMASK_spectral_Driver call basic_init() #ifdef PETSc - case (DAMASK_spectral_SolverBasicPETSC_label) - call BasicPETSC_init() + case (DAMASK_spectral_SolverBasicPETSc_label) + call basicPETSc_init() case (DAMASK_spectral_SolverAL_label) call AL_init() @@ -316,7 +312,6 @@ program DAMASK_spectral_Driver if (debugGeneral) write(6,'(a)') 'Header of result file written out' endif - !-------------------------------------------------------------------------------------------------- ! loopping over loadcases loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) @@ -335,12 +330,12 @@ program DAMASK_spectral_Driver !-------------------------------------------------------------------------------------------------- ! forwarding time timeinc_old = 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 + 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 - if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale + if (currentLoadCase == 1_pInt) then ! 1st currentLoadCase of logarithmic scale if (inc == 1_pInt) then ! 1st inc of 1st currentLoadCase of logarithmic scale - timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd + timeinc = loadCases(1)%time*(2.0_pReal**real( 1_pInt-loadCases(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd else ! not-1st inc of 1st currentLoadCase of logarithmic scale timeinc = loadCases(1)%time*(2.0_pReal**real(inc-1_pInt-loadCases(1)%incs ,pReal)) endif @@ -351,63 +346,92 @@ program DAMASK_spectral_Driver real(loadCases(currentLoadCase)%incs ,pReal)) ) endif endif - time = time + timeinc + timeinc = timeinc / 2.0_pReal**real(cutBackLevel,pReal) if(totalIncsCounter >= restartInc) then ! do calculations (otherwise just forwarding) - + stepFraction = 0_pInt !-------------------------------------------------------------------------------------------------- ! report begin of new increment - write(6,'(a)') '##################################################################' - write(6,'(A,I5.5,A,es12.5)') 'Increment ', totalIncsCounter, ' Time ',time - - select case(myspectralsolver) - - case (DAMASK_spectral_SolverBasic_label) - solres = basic_solution (& - guessmode,timeinc,timeinc_old, & - P_BC = loadCases(currentLoadCase)%P, & - F_BC = loadCases(currentLoadCase)%deformation, & - temperature_bc = loadCases(currentLoadCase)%temperature, & - rotation_BC = loadCases(currentLoadCase)%rotation) + subIncLooping: do while (stepFraction/2_pInt**cutBackLevel <1_pInt) + time = time + timeinc + stepFraction = stepFraction + 1_pInt + write(6,'(1/,a)') '###########################################################################' + write(6,'(a,es12.5'//& + ',a,'//IO_intOut(inc)//',a,'//IO_intOut(loadCases(currentLoadCase)%incs)//& + ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(2_pInt**cutBackLevel)//& + ',a,'//IO_intOut(currentLoadCase)//',a,'//IO_intOut(size(loadCases))//')') & + 'Time', time, & + 's: Increment ', inc, '/', loadCases(currentLoadCase)%incs,& + '-', stepFraction, '/', 2_pInt**cutBackLevel,& + ' of load case ', currentLoadCase,'/',size(loadCases) + write(incInfo,'(a,'//IO_intOut(totalIncsCounter)//',a,'//IO_intOut(sum(loadCases(:)%incs))//& + ',a,'//IO_intOut(stepFraction)//',a,'//IO_intOut(2_pInt**cutBackLevel)//')') & + 'Inc. ',totalIncsCounter,'/',sum(loadCases(:)%incs),& + '-',stepFraction, '/', 2_pInt**cutBackLevel + select case(myspectralsolver) + + case (DAMASK_spectral_SolverBasic_label) + solres = basic_solution (& + incInfo, guessmode,timeinc,timeinc_old, & + P_BC = loadCases(currentLoadCase)%P, & + F_BC = loadCases(currentLoadCase)%deformation, & + temperature_bc = loadCases(currentLoadCase)%temperature, & + rotation_BC = loadCases(currentLoadCase)%rotation) #ifdef PETSc - case (DAMASK_spectral_SolverBasicPETSC_label) - solres = BasicPETSC_solution (& - guessmode,timeinc,timeinc_old, & - P_BC = loadCases(currentLoadCase)%P, & - F_BC = loadCases(currentLoadCase)%deformation, & - temperature_bc = loadCases(currentLoadCase)%temperature, & - rotation_BC = loadCases(currentLoadCase)%rotation) - - case (DAMASK_spectral_SolverAL_label) - solres = AL_solution (& - guessmode,timeinc,timeinc_old, & - P_BC = loadCases(currentLoadCase)%P, & - F_BC = loadCases(currentLoadCase)%deformation, & - temperature_bc = loadCases(currentLoadCase)%temperature, & - rotation_BC = loadCases(currentLoadCase)%rotation) + case (DAMASK_spectral_SolverBasicPETSC_label) + solres = BasicPETSC_solution (& + incInfo, guessmode,timeinc,timeinc_old, & + P_BC = loadCases(currentLoadCase)%P, & + F_BC = loadCases(currentLoadCase)%deformation, & + temperature_bc = loadCases(currentLoadCase)%temperature, & + rotation_BC = loadCases(currentLoadCase)%rotation) + + case (DAMASK_spectral_SolverAL_label) + solres = AL_solution (& + incInfo, guessmode,timeinc,timeinc_old, & + P_BC = loadCases(currentLoadCase)%P, & + F_BC = loadCases(currentLoadCase)%deformation, & + temperature_bc = loadCases(currentLoadCase)%temperature, & + rotation_BC = loadCases(currentLoadCase)%rotation) #endif - end select - - write(6,'(a)') '' - write(6,'(a)') '==================================================================' + end select + cutBack = .False. + if(solres%termIll .or. .not. solres%converged) then ! no solution found + if (cutBackLevel < maxCutBack) then ! do cut back + write(6,'(/,a)') 'cut back detected' + cutBack = .True. + stepFraction = (stepFraction - 1_pInt) * 2_pInt**cutBackLevel + cutBackLevel = cutBackLevel + 1_pInt + time = time - timeinc + timeinc_old = timeinc + timeinc = timeinc/2.0_pReal + elseif (solres%termIll) then ! material point model cannot find a solution + call IO_error(850_pInt) + else + guessmode=1.0_pReal ! start guessing after first accepted (not converged) (sub)inc + endif + else + guessmode = 1.0_pReal ! start guessing after first converged (sub)inc + endif + enddo subIncLooping + cutBackLevel = max(0_pInt, cutBackLevel - 1_pInt) ! try half subincs next inc if(solres%converged) then convergedCounter = convergedCounter + 1_pInt write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') & - 'increment', totalIncsCounter, 'converged' + 'increment ', totalIncsCounter, ' converged' else write(6,'(A,'//IO_intOut(totalIncsCounter)//',A)') & - 'increment', totalIncsCounter, 'NOT converged' + 'increment ', totalIncsCounter, ' NOT converged' notConvergedCounter = notConvergedCounter + 1_pInt endif - if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency - write(6,'(a)') '' - write(6,'(a)') '... writing results to file ......................................' - write(538) materialpoint_results ! write result to file + if (mod(inc,loadCases(currentLoadCase)%outputFrequency) == 0_pInt) then ! at output frequency + write(6,'(1/,a)') '... writing results to file ......................................' + write(538) materialpoint_results ! write result to file endif - - endif ! end calculation/forwarding - guessmode = 1.0_pReal ! keep guessing along former trajectory during same currentLoadCase + else !just time forwarding + time = time + timeinc + endif ! end calculation/forwarding enddo incLooping enddo loadCaseLooping diff --git a/code/DAMASK_spectral_SolverAL.f90 b/code/DAMASK_spectral_SolverAL.f90 index a1844dbf9..cc73d3ef5 100644 --- a/code/DAMASK_spectral_SolverAL.f90 +++ b/code/DAMASK_spectral_SolverAL.f90 @@ -15,7 +15,7 @@ module DAMASK_spectral_SolverAL math_I3 use DAMASK_spectral_Utilities, only: & - solutionState + tSolutionState implicit none #include @@ -31,47 +31,48 @@ module DAMASK_spectral_SolverAL #include character (len=*), parameter, public :: & - DAMASK_spectral_SolverAL_label = 'AL' + DAMASK_spectral_SolverAL_label = 'al' !-------------------------------------------------------------------------------------------------- ! derived types - type solutionParams + type tSolutionParams real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal) :: timeinc - end type solutionParams + end type tSolutionParams - type(solutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data - DM, private :: da - SNES, private :: snes - Vec, private :: solution_vec + DM, private :: da + SNES, private :: snes + Vec, private :: solution_vec !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc - real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates - real(pReal), private, dimension(:,:,:), allocatable :: temperature + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, F_lambda_lastInc, F_lambdaDot, Fdot + real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates + real(pReal), private, dimension(:,:,:), allocatable :: temperature !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), private, dimension(3,3) :: & - F_aim = math_I3, & - F_aim_lastInc = math_I3, & - P_av - - real(pReal), private, dimension(3,3,3,3) :: & - C = 0.0_pReal, & - S = 0.0_pReal, & - C_scale = 0.0_pReal, & - S_scale = 0.0_pReal + real(pReal), private, dimension(3,3) :: & + F_aimDot, & + F_aim = math_I3, & + F_aim_lastInc = math_I3, & + P_av + character(len=1024), private :: incInfo + real(pReal), private, dimension(3,3,3,3) :: & + C = 0.0_pReal, C_lastInc = 0.0_pReal, & + S = 0.0_pReal, & + C_scale = 0.0_pReal, & + S_scale = 0.0_pReal - real(pReal), private :: err_stress, err_f, err_p - logical, private :: ForwardData - real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal + real(pReal), private :: err_stress, err_f, err_p + logical, private :: ForwardData + real(pReal), private, dimension(3,3) :: mask_stress = 0.0_pReal - contains +contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields and fills them with data, potentially from restart info @@ -93,7 +94,7 @@ subroutine AL_init() Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - debugrestart + debugRestart use numerics, only: & petsc_options @@ -118,14 +119,16 @@ subroutine AL_init() call Utilities_init() - write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverAL init -+>>>' write(6,'(a)') ' $Id: DAMASK_spectral_SolverAL.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' #include "compilation_info.f90" write(6,'(a)') '' allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) + ! allocate (Fdot,source = F_lastInc) somethin like that should be possible allocate (F_lambda_lastInc(3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (F_lambdaDot(3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) @@ -214,7 +217,8 @@ subroutine AL_init() !-------------------------------------------------------------------------------------------------- !> @brief solution for the AL scheme with internal iterations !-------------------------------------------------------------------------------------------------- - type(solutionState) function AL_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) + type(tSolutionState) function & + AL_solution(incInfoIn,guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) use numerics, only: & update_gamma @@ -227,24 +231,27 @@ subroutine AL_init() deformed_fft use IO, only: & IO_write_JobBinaryFile - use DAMASK_spectral_Utilities, only: & - boundaryCondition, & + tBoundaryCondition, & Utilities_forwardField, & + Utilities_calculateRate, & Utilities_maskedCompliance, & - Utilities_updateGamma + Utilities_updateGamma, & + cutBack use FEsolving, only: & - restartWrite + restartWrite, & + terminallyIll implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode - type(boundaryCondition), intent(in) :: P_BC,F_BC + type(tBoundaryCondition), intent(in) :: P_BC,F_BC + character(len=*), intent(in) :: incInfoIn real(pReal), dimension(3,3), intent(in) :: rotation_BC - real(pReal), dimension(3,3) :: deltaF_aim, & - F_aim_lab + real(pReal), dimension(3,3) ,save :: F_aimDot + real(pReal), dimension(3,3) :: F_aim_lab !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. real(pReal), dimension(3,3) :: temp33_Real @@ -257,6 +264,7 @@ subroutine AL_init() !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver + incInfo = incInfoIn if (restartWrite) then write(6,'(a)') 'writing converged results for restart' call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) @@ -266,31 +274,56 @@ subroutine AL_init() write (777,rec=1) C close(777) endif - AL_solution%converged =.false. + AL_solution%converged =.false. + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc) + F => xx_psc(0:8,:,:,:) + F_lambda => xx_psc(9:17,:,:,:) + + +if ( cutBack) then + F_aim = F_aim_lastInc + F_lambda= reshape(F_lambda_lastInc,[9,res(1),res(2),res(3)]) + F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + C = C_lastInc +else + !-------------------------------------------------------------------------------------------------- -! winding forward of deformation aim in loadcase system - if (F_BC%myType=='l') then ! calculate deltaF_aim from given L and current F - deltaF_aim = timeinc * F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! deltaF_aim = fDot *timeinc where applicable - deltaF_aim = timeinc * F_BC%maskFloat * F_BC%values + C_lastInc = C +!-------------------------------------------------------------------------------------------------- +! 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 + f_aimDot = F_BC%maskFloat * F_BC%values endif - temp33_Real = F_aim - F_aim = F_aim & - + guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)*timeinc/timeinc_old & - + deltaF_aim - F_aim_lastInc = temp33_Real - F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame - + f_aimDot = f_aimDot & + + guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc,timeinc_old,guessmode,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) + F_lambdaDot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc,timeinc_old,guessmode,F_lambda_lastInc,reshape(F_lambda,[3,3,res(1),res(2),res(3)])) + + F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) + F_lambda_lastInc = reshape(F_lambda,[3,3,res(1),res(2),res(3)]) + call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & + 1.0_pReal,F_lastInc,coordinates) + endif + F_aim = F_aim + f_aimDot * timeinc + !-------------------------------------------------------------------------------------------------- ! update local deformation gradient and coordinates - deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc) - F => xx_psc(0:8,:,:,:) - F_lambda => xx_psc(9:17,:,:,:) - call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F) - call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lambda_lastInc,F_lambda) + ! deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) + + + F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)]) + F_lambda = reshape(Utilities_forwardField(timeinc,F_aim,F_lambda_lastInc,F_lambdadot),[9,res(1),res(2),res(3)]) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc) - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) @@ -305,6 +338,11 @@ subroutine AL_init() call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr_psc) call SNESGetConvergedReason(snes,reason,ierr_psc) + + AL_solution%termIll = terminallyIll + terminallyIll = .false. + + if (reason > 0 ) AL_solution%converged = .true. end function AL_solution @@ -326,15 +364,18 @@ subroutine AL_init() wgt use DAMASK_spectral_Utilities, only: & field_real, & - Utilities_forwardFFT, & + Utilities_FFTforward, & Utilities_fourierConvolution, & - Utilities_backwardFFT, & + Utilities_FFTbackward, & Utilities_constitutiveResponse + use IO, only: IO_intOut implicit none integer(pInt) :: i,j,k - real(pReal), dimension(3,3) :: temp33_Real + integer(pInt), save :: callNo = 3_pInt, reportIter = 0_pInt + real(pReal), dimension(3,3) :: temp33_Real + logical :: report DMDALocalInfo :: in(DMDA_LOCAL_INFO_SIZE) PetscScalar, target :: x_scal(3,3,2,XG_RANGE,YG_RANGE,ZG_RANGE) @@ -355,14 +396,21 @@ subroutine AL_init() !-------------------------------------------------------------------------------------------------- ! report begin of new iteration - write(6,'(a)') '' - write(6,'(a)') '==================================================================' - write(6,'(4(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax, ' | # Func. calls = ',nfuncs + if (iter == 0 .and. callNo>2) then + callNo = 0_pInt + reportIter = 0_pInt + endif + if (callNo == 0 .or. mod(callNo,2) == 1_pInt) then + write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iter. ', itmin, '<',reportIter, '≤', itmax write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& math_transpose33(F_aim) - + reportIter = reportIter + 1_pInt + endif + callNo = callNo +1_pInt !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response + call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, & residual_F,C,P_av,ForwardData,params%rotation_BC) ForwardData = .False. @@ -384,9 +432,9 @@ subroutine AL_init() !-------------------------------------------------------------------------------------------------- ! doing Fourier transform - call Utilities_forwardFFT() + call Utilities_FFTforward() call Utilities_fourierConvolution(math_rotate_backward33(F_aim,params%rotation_BC)) - call Utilities_backwardFFT() + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual @@ -398,7 +446,7 @@ subroutine AL_init() ! calculating errors err_f = wgt*sqrt(sum(residual_F_lambda**2.0_pReal)) err_p = wgt*sqrt(sum((residual_F - residual_F_lambda)**2.0_pReal)) - + end subroutine AL_formResidual !-------------------------------------------------------------------------------------------------- @@ -440,8 +488,7 @@ subroutine AL_init() write(6,'(a,es14.7)') 'error stress BC = ', err_stress/min(maxval(abs(P_av))*err_stress_tolrel,err_stress_tolabs) write(6,'(a,es14.7)') 'error F = ', err_f/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_f_tol write(6,'(a,es14.7)') 'error P = ', err_p/sqrt(sum((F_aim-math_I3)*(F_aim-math_I3)))/err_p_tol - return - + write(6,'(/,a)') '==========================================================================' end subroutine AL_converged !-------------------------------------------------------------------------------------------------- diff --git a/code/DAMASK_spectral_SolverBasic.f90 b/code/DAMASK_spectral_SolverBasic.f90 index a3c2ad7ea..c7360f923 100644 --- a/code/DAMASK_spectral_SolverBasic.f90 +++ b/code/DAMASK_spectral_SolverBasic.f90 @@ -17,7 +17,7 @@ module DAMASK_spectral_SolverBasic math_I3 use DAMASK_spectral_Utilities, only: & - solutionState + tSolutionState implicit none character (len=*), parameter, public :: & @@ -25,7 +25,7 @@ module DAMASK_spectral_SolverBasic !-------------------------------------------------------------------------------------------------- ! pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, P + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, Fdot, P real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal), private, dimension(:,:,:), allocatable :: temperature @@ -35,7 +35,7 @@ module DAMASK_spectral_SolverBasic F_aim = math_I3, & F_aim_lastInc = math_I3 real(pReal), private,dimension(3,3,3,3) :: & - C = 0.0_pReal + C = 0.0_pReal, C_lastInc = 0.0_pReal contains @@ -60,28 +60,26 @@ subroutine basic_init() Utilities_init, & Utilities_constitutiveResponse, & Utilities_updateGamma, & - debugrestart + debugRestart use mesh, only: & res, & geomdim implicit none - integer(pInt) :: i,j,k real(pReal), dimension(3,3) :: temp33_Real call Utilities_Init() - write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" write(6,'(a)') '' - allocate (F ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (F_lastInc ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) + allocate (Fdot ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (P ( 3,3,res(1), res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) @@ -139,7 +137,8 @@ end subroutine basic_init !-------------------------------------------------------------------------------------------------- !> @brief solution for the basic scheme with internal iterations !-------------------------------------------------------------------------------------------------- -type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) +type(tSolutionState) function & + basic_solution(incInfo,guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) use numerics, only: & itmax, & @@ -153,39 +152,40 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F use mesh, only: & res,& geomdim, & - deformed_fft + deformed_fft, & + wgt use IO, only: & IO_write_JobBinaryFile, & IO_intOut use DAMASK_spectral_Utilities, only: & - boundaryCondition, & + tBoundaryCondition, & field_real, & Utilities_forwardField, & Utilities_maskedCompliance, & - Utilities_forwardFFT, & + Utilities_FFTforward, & Utilities_divergenceRMS, & Utilities_fourierConvolution, & - Utilities_backwardFFT, & + Utilities_FFTbackward, & Utilities_updateGamma, & - Utilities_constitutiveResponse + Utilities_constitutiveResponse, & + Utilities_calculateRate use FEsolving, only: & restartWrite, & terminallyIll - + use DAMASK_spectral_Utilities, only: cutBack implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode - type(boundaryCondition), intent(in) :: P_BC,F_BC + type(tBoundaryCondition), intent(in) :: P_BC,F_BC + character(len=*), intent(in) :: incInfo real(pReal), dimension(3,3), intent(in) :: rotation_BC - - real(pReal), dimension(3,3,3,3) :: S - real(pReal), dimension(3,3) :: deltaF_aim, & - F_aim_lab, & + real(pReal), dimension(3,3), save :: f_aimDot = 0.0_pReal + real(pReal), dimension(3,3) :: F_aim_lab, & F_aim_lab_lastIter, & P_av !-------------------------------------------------------------------------------------------------- @@ -208,69 +208,80 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F close(777) endif + + if ( cutBack) then + F_aim = F_aim_lastInc + F = F_lastInc + C = C_lastInc +else + C_lastInc = C !-------------------------------------------------------------------------------------------------- -! winding forward of deformation aim in loadcase system - if (F_BC%myType=='l') then ! calculate deltaF_aim from given L and current F - deltaF_aim = timeinc * F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! deltaF_aim = fDot *timeinc where applicable - deltaF_aim = timeinc * F_BC%maskFloat * F_BC%values +! 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 + f_aimDot = F_BC%maskFloat * F_BC%values + endif + f_aimDot = f_aimDot & + + guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + +!-------------------------------------------------------------------------------------------------- +! update coordinates and rate and forward last inc + call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & + 1.0_pReal,F_lastInc,coordinates) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc,timeinc_old,guessmode,F_lastInc,F) + F_lastInc = F endif - temp33_Real = F_aim - F_aim = F_aim & - + guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)*timeinc/timeinc_old & - + deltaF_aim - F_aim_lastInc = temp33_Real - F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame - -!-------------------------------------------------------------------------------------------------- -! update local deformation gradient and coordinates - deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) - call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F) - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) - + F_aim = F_aim + f_aimDot * timeinc + F = Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot) !I thin F aim should be rotated here +print*, 'F', sum(sum(sum(F,dim=5),dim=4),dim=3)*wgt !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) if (update_gamma) call Utilities_updateGamma(C) + iter = 0_pInt - ForwardData = .True. convergenceLoop: do while(iter < itmax) iter = iter + 1_pInt !-------------------------------------------------------------------------------------------------- ! report begin of new iteration - write(6,'(a)') '' - write(6,'(a)') '==================================================================' - write(6,'(3(a,'//IO_intOut(itmax)//'))') ' Iter.', itmin, '<',iter, '<', itmax + 1_pInt + write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iter. ', itmin, '<',iter, '≤', itmax write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =', & math_transpose33(F_aim) F_aim_lab_lastIter = math_rotate_backward33(F_aim,rotation_BC) !-------------------------------------------------------------------------------------------------- ! evaluate constitutive response + basic_solution%termIll = .false. call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& P,C,P_av,ForwardData,rotation_BC) basic_solution%termIll = terminallyIll + terminallyIll = .False. ForwardData = .False. !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC%values))) !S = 0.0 for no bc - err_stress = maxval(abs(P_BC%maskFloat * (P_av - P_BC%values))) ! mask = 0.0 for no bc - F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame + F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC%values))) ! S = 0.0 for no bc + err_stress = maxval(abs(P_BC%maskFloat * (P_av - P_BC%values))) ! mask = 0.0 for no bc + F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame !-------------------------------------------------------------------------------------------------- ! updated deformation gradient using fix point algorithm of basic scheme field_real = 0.0_pReal field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(P,[res(1),res(2),res(3),3,3],& order=[4,5,1,2,3]) ! field real has a different order - call Utilities_forwardFFT() + call Utilities_FFTforward() err_div = Utilities_divergenceRMS() call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) - call Utilities_backwardFFT() + call Utilities_FFTbackward() F = F - reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization basic_solution%converged = basic_Converged(err_div,P_av,err_stress,P_av) + write(6,'(/,a)') '==========================================================================' if ((basic_solution%converged .and. iter > itmin) .or. basic_solution%termIll) exit enddo convergenceLoop @@ -320,6 +331,7 @@ logical function basic_Converged(err_div,pAvgDiv,err_stress,pAvgStress) end function basic_Converged + subroutine basic_destroy() use DAMASK_spectral_Utilities, only: & @@ -330,6 +342,9 @@ subroutine basic_destroy() end subroutine basic_destroy + + + end module DAMASK_spectral_SolverBasic diff --git a/code/DAMASK_spectral_SolverBasicPETSC.f90 b/code/DAMASK_spectral_SolverBasicPETSC.f90 index b2af3555c..b90e71035 100644 --- a/code/DAMASK_spectral_SolverBasicPETSC.f90 +++ b/code/DAMASK_spectral_SolverBasicPETSC.f90 @@ -15,8 +15,8 @@ module DAMASK_spectral_SolverBasicPETSC math_I3 use DAMASK_spectral_Utilities, only: & - solutionState - + tSolutionState + implicit none #include #include @@ -35,12 +35,12 @@ module DAMASK_spectral_SolverBasicPETSC !-------------------------------------------------------------------------------------------------- ! derived types - type solutionParams + type tSolutionParams real(pReal), dimension(3,3) :: P_BC, rotation_BC real(pReal) :: timeinc - end type solutionParams + end type tSolutionParams - type(solutionParams), private :: params + type(tSolutionParams), private :: params !-------------------------------------------------------------------------------------------------- ! PETSc data @@ -50,7 +50,7 @@ module DAMASK_spectral_SolverBasicPETSC !-------------------------------------------------------------------------------------------------- ! common pointwise data - real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F_lastInc, Fdot real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates real(pReal), private, dimension(:,:,:), allocatable :: temperature @@ -60,9 +60,9 @@ module DAMASK_spectral_SolverBasicPETSC F_aim = math_I3, & F_aim_lastInc = math_I3, & P_av - + character(len=1024), private :: incInfo real(pReal), private, dimension(3,3,3,3) :: & - C = 0.0_pReal, & + C = 0.0_pReal, C_lastInc= 0.0_pReal, & S = 0.0_pReal real(pReal), private :: err_stress, err_div @@ -78,52 +78,51 @@ subroutine BasicPETSC_init() use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use IO, only: & - IO_read_JobBinaryFile, & - IO_write_JobBinaryFile + use IO, only: & + IO_read_JobBinaryFile, & + IO_write_JobBinaryFile - use FEsolving, only: & - restartInc + use FEsolving, only: & + restartInc - use DAMASK_interface, only: & - getSolverJobName + use DAMASK_interface, only: & + getSolverJobName - use DAMASK_spectral_Utilities, only: & - Utilities_init, & - Utilities_constitutiveResponse, & - Utilities_updateGamma, & - debugrestart + use DAMASK_spectral_Utilities, only: & + Utilities_init, & + Utilities_constitutiveResponse, & + Utilities_updateGamma, & + debugRestart - use numerics, only: & - petsc_options + use numerics, only: & + petsc_options - use mesh, only: & - res, & - geomdim, & - mesh_NcpElems + use mesh, only: & + res, & + geomdim, & + mesh_NcpElems - use math, only: & - math_invSym3333 + use math, only: & + math_invSym3333 - implicit none + implicit none + integer(pInt) :: i,j,k + real(pReal), dimension(:,:,:,:,:), allocatable :: P - integer(pInt) :: i,j,k - real(pReal), dimension(:,:,:,:,:), allocatable :: P - - PetscErrorCode :: ierr_psc - PetscObject :: dummy - PetscMPIInt :: rank - PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:) + PetscErrorCode :: ierr_psc + PetscObject :: dummy + PetscMPIInt :: rank + PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:) call Utilities_init() - write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectral_solverBasicPETSC init -+>>>' + write(6,'(/,a)') ' <<<+- DAMASK_spectral_solverBasicPETSC init -+>>>' write(6,'(a)') ' $Id: DAMASK_spectral_SolverBasicPETSC.f90 1654 2012-08-03 09:25:48Z MPIE\m.diehl $' #include "compilation_info.f90" write(6,'(a)') '' allocate (F_lastInc (3,3, res(1), res(2),res(3)), source = 0.0_pReal) + allocate (Fdot (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (P (3,3, res(1), res(2),res(3)), source = 0.0_pReal) allocate (coordinates( res(1), res(2),res(3),3), source = 0.0_pReal) allocate (temperature( res(1), res(2),res(3)), source = 0.0_pReal) @@ -199,8 +198,8 @@ subroutine BasicPETSC_init() !-------------------------------------------------------------------------------------------------- !> @brief solution for the Basic PETSC scheme with internal iterations !-------------------------------------------------------------------------------------------------- - type(solutionState) function BasicPETSC_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) - + type(tSolutionState) function & + basicPETSc_solution(incInfoIn,guessmode,timeinc,timeinc_old,P_BC,F_BC,temperature_bc,rotation_BC) use numerics, only: & update_gamma use math, only: & @@ -212,23 +211,25 @@ subroutine BasicPETSC_init() deformed_fft use IO, only: & IO_write_JobBinaryFile - use DAMASK_spectral_Utilities, only: & - boundaryCondition, & + tBoundaryCondition, & + Utilities_calculateRate, & Utilities_forwardField, & Utilities_maskedCompliance, & - Utilities_updateGamma - + Utilities_updateGamma, & + cutBack use FEsolving, only: & - restartWrite - + restartWrite, & + terminallyIll implicit none !-------------------------------------------------------------------------------------------------- ! input data for solution real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode - type(boundaryCondition), intent(in) :: P_BC,F_BC + type(tBoundaryCondition), intent(in) :: P_BC,F_BC real(pReal), dimension(3,3), intent(in) :: rotation_BC - real(pReal), dimension(3,3) :: deltaF_aim, & + character(len=*), intent(in) :: incInfoIn + real(pReal), dimension(3,3),save :: F_aimDot=0.0_pReal + real(pReal), dimension(3,3) :: & F_aim_lab !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. @@ -237,11 +238,12 @@ subroutine BasicPETSC_init() !-------------------------------------------------------------------------------------------------- ! PetscScalar, pointer :: xx_psc(:,:,:,:), F(:,:,:,:) - PetscErrorCode ierr_psc - SNESConvergedReason reason + PetscErrorCode :: ierr_psc + SNESConvergedReason :: reason !-------------------------------------------------------------------------------------------------- ! restart information for spectral solver + incInfo = incInfoIn if (restartWrite) then write(6,'(a)') 'writing converged results for restart' call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) @@ -251,29 +253,43 @@ subroutine BasicPETSC_init() write (777,rec=1) C close(777) endif - BasicPETSC_solution%converged =.false. + call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc) + F => xx_psc(0:8,:,:,:) + +if ( cutBack) then + F_aim = F_aim_lastInc + + F = reshape(F_lastInc,[9,res(1),res(2),res(3)]) + C = C_lastInc +else + + C_lastInc = C + !-------------------------------------------------------------------------------------------------- -! winding forward of deformation aim in loadcase system - if (F_BC%myType=='l') then ! calculate deltaF_aim from given L and current F - deltaF_aim = timeinc * F_BC%maskFloat * math_mul33x33(F_BC%values, F_aim) - elseif(F_BC%myType=='fdot') then ! deltaF_aim = fDot *timeinc where applicable - deltaF_aim = timeinc * F_BC%maskFloat * F_BC%values +! 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 + f_aimDot = F_BC%maskFloat * F_BC%values endif - temp33_Real = F_aim - F_aim = F_aim & - + guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)*timeinc/timeinc_old & - + deltaF_aim - F_aim_lastInc = temp33_Real - F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame - + f_aimDot = f_aimDot & + + guessmode * P_BC%maskFloat * (F_aim - F_aim_lastInc)/timeinc_old + F_aim_lastInc = F_aim + !-------------------------------------------------------------------------------------------------- -! update local deformation gradient and coordinates - deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) - call DMDAVecGetArrayF90(da,solution_vec,xx_psc,ierr_psc) - F => xx_psc(0:8,:,:,:) - call Utilities_forwardField(deltaF_aim,timeinc,timeinc_old,guessmode,F_lastInc,F) - call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc) - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) +! update coordinates and rate and forward last inc + call deformed_fft(res,geomdim,math_rotate_backward33(F_aim_lastInc,rotation_BC), & + 1.0_pReal,F_lastInc,coordinates) + Fdot = Utilities_calculateRate(math_rotate_backward33(f_aimDot,rotation_BC), & + timeinc,timeinc_old,guessmode,F_lastInc,reshape(F,[3,3,res(1),res(2),res(3)])) + F_lastInc = reshape(F,[3,3,res(1),res(2),res(3)]) + endif + F_aim = F_aim + f_aimDot * timeinc + + + F = reshape(Utilities_forwardField(timeinc,F_aim,F_lastInc,Fdot),[9,res(1),res(2),res(3)]) + call DMDAVecRestoreArrayF90(da,solution_vec,xx_psc,ierr_psc) + call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),1.0_pReal,F_lastInc,coordinates) !-------------------------------------------------------------------------------------------------- ! update stiffness (and gamma operator) @@ -288,6 +304,9 @@ subroutine BasicPETSC_init() call SNESSolve(snes,PETSC_NULL_OBJECT,solution_vec,ierr_psc) call SNESGetConvergedReason(snes,reason,ierr_psc) + basicPETSc_solution%termIll = terminallyIll + terminallyIll = .false. + BasicPETSC_solution%converged =.false. if (reason > 0 ) BasicPETSC_solution%converged = .true. end function BasicPETSC_solution @@ -309,12 +328,12 @@ subroutine BasicPETSC_init() wgt use DAMASK_spectral_Utilities, only: & field_real, & - Utilities_forwardFFT, & + Utilities_FFTforward, & + Utilities_FFTbackward, & Utilities_fourierConvolution, & - Utilities_backwardFFT, & Utilities_constitutiveResponse, & Utilities_divergenceRMS - + use IO, only : IO_intOut implicit none real(pReal), dimension(3,3) :: F_aim_lab_lastIter, F_aim_lab @@ -327,7 +346,6 @@ subroutine BasicPETSC_init() PetscInt :: iter, nfuncs PetscObject :: dummy PetscErrorCode :: ierr_psc - F => x_scal(:,:,:,:,:) residual_F => f_scal(:,:,:,:,:) @@ -336,9 +354,8 @@ subroutine BasicPETSC_init() !-------------------------------------------------------------------------------------------------- ! report begin of new iteration - write(6,'(a)') '' - write(6,'(a)') '==================================================================' - write(6,'(4(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax, ' | # Func. calls = ',nfuncs + write(6,'(/,a,3(a,'//IO_intOut(itmax)//'))') trim(incInfo), & + ' @ Iter. ', itmin, '<',iter, '≤', itmax write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& math_transpose33(F_aim) F_aim_lab_lastIter = math_rotate_backward33(F_aim,params%rotation_BC) @@ -347,7 +364,7 @@ subroutine BasicPETSC_init() ! evaluate constitutive response call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,params%timeinc, & residual_F,C,P_av,ForwardData,params%rotation_BC) - ForwardData = .False. + ForwardData = .false. !-------------------------------------------------------------------------------------------------- ! stress BC handling @@ -360,15 +377,15 @@ subroutine BasicPETSC_init() field_real = 0.0_pReal field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = reshape(residual_F,[res(1),res(2),res(3),3,3],& order=[4,5,1,2,3]) ! field real has a different order - call Utilities_forwardFFT() + call Utilities_FFTforward() err_div = Utilities_divergenceRMS() call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) - call Utilities_backwardFFT() + call Utilities_FFTbackward() !-------------------------------------------------------------------------------------------------- ! constructing residual residual_F = reshape(field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),shape(F),order=[3,4,5,1,2]) - + write(6,'(/,a)') '==========================================================================' end subroutine BasicPETSC_formResidual !-------------------------------------------------------------------------------------------------- diff --git a/code/DAMASK_spectral_Utilities.f90 b/code/DAMASK_spectral_Utilities.f90 index 89d28bad2..4f2fe5edb 100644 --- a/code/DAMASK_spectral_Utilities.f90 +++ b/code/DAMASK_spectral_Utilities.f90 @@ -33,7 +33,7 @@ module DAMASK_spectral_Utilities real(pReal), private, dimension(3,3,3,3) :: C_ref real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real - + logical, public :: cutBack =.false. !-------------------------------------------------------------------------------------------------- ! debug fftw type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back @@ -45,7 +45,7 @@ module DAMASK_spectral_Utilities type(C_PTR), private :: plan_divergence real(pReal), private, dimension(:,:,:,:), pointer :: divergence_real complex(pReal), private, dimension(:,:,:,:), pointer :: divergence_fourier - real(pReal), dimension(:,:,:,:), allocatable :: divergence_post + real(pReal), private, dimension(:,:,:,:), allocatable :: divergence_post !-------------------------------------------------------------------------------------------------- !variables controlling debugging @@ -53,18 +53,18 @@ module DAMASK_spectral_Utilities !-------------------------------------------------------------------------------------------------- ! derived types - type solutionState - logical :: converged = .false. + type tSolutionState + logical :: converged = .true. logical :: regrid = .false. logical :: termIll = .false. - end type solutionState + end type tSolutionState - type boundaryCondition + type tBoundaryCondition real(pReal), dimension(3,3) :: values = 0.0_pReal real(pReal), dimension(3,3) :: maskFloat = 0.0_pReal logical, dimension(3,3) :: maskLogical = .false. character(len=64) :: myType = 'None' - end type boundaryCondition + end type tBoundaryCondition contains @@ -105,7 +105,7 @@ subroutine Utilities_init() write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectralSolver Utilities init -+>>>' + write(6,'(a)') ' <<<+- DAMASK_Utilities init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" write(6,'(a)') '' @@ -238,7 +238,7 @@ end subroutine Utilities_updateGamma !> In case of debugging the FFT, also one component of the tensor (specified by row and column) !> is independetly transformed complex to complex and compared to the whole tensor transform !-------------------------------------------------------------------------------------------------- -subroutine Utilities_forwardFFT(row,column) +subroutine Utilities_FFTforward(row,column) use mesh, only : & virt_dim use math, only: & @@ -286,7 +286,7 @@ subroutine Utilities_forwardFFT(row,column) if(res(3)>1_pInt) & field_fourier (1:res1_red,1:res(2), res(3)/2_pInt+1_pInt,1:3,1:3)& = cmplx(0.0_pReal,0.0_pReal,pReal) -end subroutine Utilities_forwardFFT +end subroutine Utilities_FFTforward !-------------------------------------------------------------------------------------------------- @@ -296,7 +296,7 @@ end subroutine Utilities_forwardFFT !> is independetly transformed complex to complex and compared to the whole tensor transform !> results is weighted by number of points stored in wgt !-------------------------------------------------------------------------------------------------- -subroutine Utilities_backwardFFT(row,column) +subroutine Utilities_FFTbackward(row,column) implicit none integer(pInt), intent(in), optional :: row, column @@ -335,7 +335,7 @@ subroutine Utilities_backwardFFT(row,column) endif field_real = field_real * wgt -end subroutine Utilities_backwardFFT +end subroutine Utilities_FFTbackward !-------------------------------------------------------------------------------------------------- @@ -353,9 +353,8 @@ subroutine Utilities_fourierConvolution(fieldAim) integer(pInt) :: i, j, k, l, m, n, o complex(pReal), dimension(3,3) :: temp33_complex - write(6,'(a)') '' - write(6,'(a)') '... doing convolution .................' - write(6,'(a)') '' + + write(6,'(/,a)') '... doing convolution .....................................................' !-------------------------------------------------------------------------------------------------- ! to the actual spectral method calculation (mechanical equilibrium) @@ -398,8 +397,7 @@ real(pReal) function Utilities_divergenceRMS() err_div_max, err_real_div_max complex(pReal), dimension(3) :: temp3_complex - write(6,'(a)') '' - write(6,'(a)') '... calculating divergence .................' + write(6,'(/,a)') '... calculating divergence ................................................' !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space @@ -477,7 +475,7 @@ function Utilities_maskedCompliance(rot_BC,mask_stressVector,C) allocate (s_reduced(size_reduced,size_reduced), source =0.0_pReal) allocate (sTimesC(size_reduced,size_reduced), source =0.0_pReal) - C_lastInc = math_rotate_forward3333(C*0.95_pReal+0.5_pReal*C_ref,rot_BC) ! calculate stiffness from former inc + C_lastInc = math_rotate_forward3333(C*0.75_pReal+0.05_pReal*C_ref,rot_BC) ! calculate stiffness from former inc temp99_Real = math_Plain3333to99(C_lastInc) k = 0_pInt ! build reduced stiffness do n = 1_pInt,9_pInt @@ -537,77 +535,99 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti implicit none real(pReal), dimension(res(1),res(2),res(3)) :: temperature - real(pReal), dimension(res(1),res(2),res(3),3) :: coordinates + real(pReal), dimension(res(1),res(2),res(3),3), intent(in) :: coordinates - real(pReal), dimension(3,3,res(1),res(2),res(3)) :: F,F_lastInc, P + real(pReal), dimension(3,3,res(1),res(2),res(3)), intent(in) :: F,F_lastInc + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: P real(pReal) :: timeinc - logical :: ForwardData + logical, intent(in) :: forwardData integer(pInt) :: i, j, k, ielem - integer(pInt) :: CPFEM_mode + integer(pInt) :: calcMode, collectMode real(pReal), dimension(3,3,3,3) :: dPdF, C real(pReal), dimension(6) :: sigma ! cauchy stress real(pReal), dimension(6,6) :: dsde - real(pReal), dimension(3,3) :: P_av, rotation_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC + real(pReal), dimension(3,3) :: P_av - write(6,'(a)') '' - write(6,'(a)') '... evaluating constitutive response .................' - write(6,'(a)') '' - - if (ForwardData) then - CPFEM_mode = 1_pInt + write(6,'(/,a,/)') '... evaluating constitutive response ......................................' + if (forwardData) then + calcMode = 1_pInt + collectMode = 4_pInt else - CPFEM_mode = 2_pInt + calcMode = 2_pInt + collectMode = 3_pInt endif - + if (cutBack) then + calcMode = 2_pInt + collectMode = 5_pInt + endif + print*, 'collect mode', collectMode + print*, 'calc mode', calcMode ielem = 0_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt - call CPFEM_general(3_pInt,& ! collect cycle + call CPFEM_general(collectMode,& ! collect cycle coordinates(i,j,k,1:3), F_lastInc(1:3,1:3,i,j,k),F(1:3,1:3,i,j,k), & temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) + collectMode = 3_pInt enddo; enddo; enddo - P = 0.0_pReal ! needed because of the padding for FFTW + P = 0.0_pReal ! needed because of the padding for FFTW C = 0.0_pReal ielem = 0_pInt call debug_reset() do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) ielem = ielem + 1_pInt - call CPFEM_general(CPFEM_mode,& ! first element in first iteration retains CPFEM_mode 1, - coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) + call CPFEM_general(calcMode,& ! first element in first iteration retains CPFEM_mode 1, + coordinates(i,j,k,1:3),F_lastInc(1:3,1:3,i,j,k), F(1:3,1:3,i,j,k), & ! others get 2 (saves winding forward effort) temperature(i,j,k),timeinc,ielem,1_pInt,sigma,dsde,P(1:3,1:3,i,j,k),dPdF) - CPFEM_mode = 2_pInt + calcMode = 2_pInt C = C + dPdF enddo; enddo; enddo call debug_info() P_av = math_rotate_forward33(sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt,rotation_BC) !average of P rotated restartWrite = .false. - + cutBack = .false. - write (6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') ' Piola-Kirchhoff stress / MPa =',& + write(6,'(a,/,3(3(2x,f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& math_transpose33(P_av)/1.e6_pReal C = C * wgt end subroutine Utilities_constitutiveResponse -subroutine Utilities_forwardField(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field) +function Utilities_calculateRate(delta_aim,timeinc,timeinc_old,guessmode,field_lastInc,field) implicit none real(pReal), intent(in), dimension(3,3) :: delta_aim real(pReal), intent(in) :: timeinc, timeinc_old, guessmode - real(pReal), intent(inout), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field + real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,field + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: Utilities_calculateRate if (guessmode == 1.0_pReal) then - field = field + (field-field_lastInc) * timeinc/timeinc_old - field_lastInc = (field + field_lastInc * timeinc/timeinc_old) /(1.0_pReal + timeinc/timeinc_old) + Utilities_calculateRate = (field-field_lastInc) / timeinc_old else - field_lastInc = field - field = field + spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3)) + Utilities_calculateRate = spread(spread(spread(delta_aim,3,res(1)),4,res(2)),5,res(3))/timeinc endif -end subroutine Utilities_forwardField +end function Utilities_calculateRate + + +function Utilities_forwardField(timeinc,aim,field_lastInc,rate) + implicit none + real(pReal), intent(in) :: timeinc + real(pReal), intent(in), dimension(3,3) :: aim + real(pReal), intent(in), dimension(3,3,res(1),res(2),res(3)) :: field_lastInc,rate + real(pReal), dimension(3,3,res(1),res(2),res(3)) :: Utilities_forwardField + real(pReal), dimension(3,3) :: fieldDiff + Utilities_forwardField = field_lastInc + rate*timeinc + fieldDiff = sum(sum(sum(Utilities_forwardField,dim=5),dim=4),dim=3)*wgt - aim + Utilities_forwardField = Utilities_forwardField - & + spread(spread(spread(fieldDiff,3,res(1)),4,res(2)),5,res(3)) + +end function Utilities_forwardField + real(pReal) function Utilities_getFilter(k) use numerics, only: & @@ -634,6 +654,7 @@ real(pReal) function Utilities_getFilter(k) end function Utilities_getFilter + subroutine Utilities_destroy() implicit none @@ -648,5 +669,6 @@ subroutine Utilities_destroy() call fftw_destroy_plan(plan_backward) end subroutine Utilities_destroy - + + end module DAMASK_spectral_Utilities diff --git a/code/numerics.f90 b/code/numerics.f90 index fd431ef28..50682abfa 100644 --- a/code/numerics.f90 +++ b/code/numerics.f90 @@ -25,20 +25,23 @@ module numerics use prec, only: pInt, pReal implicit none -character(len=64), parameter, private ::& +character(len=64), parameter, private :: & numerics_configFile = 'numerics.config' !< name of configuration file -integer(pInt) :: iJacoStiffness = 1_pInt, & !< frequency of stiffness update +integer(pInt), protected, public :: & + iJacoStiffness = 1_pInt, & !< frequency of stiffness update iJacoLpresiduum = 1_pInt, & !< frequency of Jacobian update of residuum in Lp nHomog = 20_pInt, & !< homogenization loop limit (only for debugging info, loop limit is determined by "subStepMinHomog") nMPstate = 10_pInt, & !< materialpoint state loop limit nCryst = 20_pInt, & !< crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nState = 10_pInt, & !< state loop limit nStress = 40_pInt, & !< stress loop limit - pert_method = 1_pInt, & !< method used in perturbation technique for tangent - numerics_integrationMode = 0_pInt !< integrationMode 1 = central solution ; integrationMode 2 = perturbation, Default 0: undefined, is not read from file -integer(pInt), dimension(2) :: numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states -real(pReal) :: relevantStrain = 1.0e-7_pReal, & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) + pert_method = 1_pInt !< method used in perturbation technique for tangent +integer(pInt) :: numerics_integrationMode = 0_pInt !< integrationMode 1 = central solution ; integrationMode 2 = perturbation, Default 0: undefined, is not read from file +integer(pInt), dimension(2) , protected, public :: & + numerics_integrator = 1_pInt !< method used for state integration (central & perturbed state), Default 1: fix-point iteration for both states +real(pReal), protected, public :: & + relevantStrain = 1.0e-7_pReal, & !< strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) defgradTolerance = 1.0e-7_pReal, & !< deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) pert_Fg = 1.0e-7_pReal, & !< strain perturbation for FEM Jacobi subStepMinCryst = 1.0e-3_pReal, & !< minimum (relative) size of sub-step allowed during cutback in crystallite @@ -65,39 +68,49 @@ real(pReal) :: relevantStrain = 1.0e-7_pReal, & maxVolDiscr_RGC = 1.0e-5_pReal, & !< threshold of maximum volume discrepancy allowed volDiscrMod_RGC = 1.0e+12_pReal, & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint) volDiscrPow_RGC = 5.0_pReal !< powerlaw penalty for volume discrepancy -logical :: analyticJaco = .false. !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations +logical, protected, public :: & + analyticJaco = .false. !< use analytic Jacobian or perturbation, Default .false.: calculate Jacobian using perturbations !* Random seeding parameters -integer(pInt) :: fixedSeed = 0_pInt !< fixed seeding for pseudo-random number generator, Default 0: use random seed +integer(pInt), protected, public :: & + fixedSeed = 0_pInt !< fixed seeding for pseudo-random number generator, Default 0: use random seed !* OpenMP variable -integer(pInt) :: DAMASK_NumThreadsInt = 0_pInt !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive +integer(pInt), protected, public :: & + DAMASK_NumThreadsInt = 0_pInt !< value stored in environment variable DAMASK_NUM_THREADS, set to zero if no OpenMP directive !* spectral parameters: #ifdef Spectral -real(pReal) :: err_div_tol = 0.1_pReal, & !< Div(P)/avg(P)*meter +real(pReal), protected, public :: & + err_div_tol = 0.1_pReal, & !< Div(P)/avg(P)*meter err_stress_tolrel = 0.01_pReal, & !< relative tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress err_stress_tolabs = huge(1.0_pReal), & !< absolute tolerance for fullfillment of stress BC, Default: 0.01 allowing deviation of 1% of maximum stress err_f_tol = 1e-6_pReal, & err_p_tol = 1e-5_pReal, & fftw_timelimit = -1.0_pReal, & !< sets the timelimit of plan creation for FFTW, see manual on www.fftw.org, Default -1.0: disable timelimit rotation_tol = 1.0e-12_pReal !< tolerance of rotation specified in loadcase, Default 1.0e-12: first guess -character(len=64) :: fftw_plan_mode = 'FFTW_PATIENT', & !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag +character(len=64), private :: & + fftw_plan_mode = 'FFTW_PATIENT' !< reads the planing-rigor flag, see manual on www.fftw.org, Default FFTW_PATIENT: use patient planner flag +character(len=64), protected, public :: & myspectralsolver = 'basic' , & !< spectral solution method myfilter = 'none' !< spectral filtering method -character(len=1024) :: petsc_options = '-snes_type ngmres -snes_ngmres_anderson -snes_view' -integer(pInt) :: fftw_planner_flag = 32_pInt, & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw +character(len=1024), protected, public :: & + petsc_options = '-snes_type ngmres -snes_ngmres_anderson -snes_view' +integer(pInt), protected, public :: & + fftw_planner_flag = 32_pInt, & !< conversion of fftw_plan_mode to integer, basically what is usually done in the include file of fftw itmax = 20_pInt, & !< maximum number of iterations itmin = 2_pInt, & !< minimum number of iterations maxCutBack = 3_pInt !< max number of cut backs -logical :: memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate +logical, protected , public :: & + memory_efficient = .true., & !< for fast execution (pre calculation of gamma_hat), Default .true.: do not precalculate divergence_correction = .false., & !< correct divergence calculation in fourier space, Default .false.: no correction update_gamma = .false. !< update gamma operator with current stiffness, Default .false.: use initial stiffness + #endif + public :: numerics_init + +contains - -CONTAINS - !******************************************* ! initialization subroutine !******************************************* @@ -256,9 +269,9 @@ subroutine numerics_init case ('fftw_timelimit') fftw_timelimit = IO_floatValue(line,positions,2_pInt) case ('fftw_plan_mode') - fftw_plan_mode = IO_stringValue(line,positions,2_pInt) + fftw_plan_mode = IO_lc(IO_stringValue(line,positions,2_pInt)) case ('myfilter') - myfilter = IO_stringValue(line,positions,2_pInt) + myfilter = IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rotation_tol') rotation_tol = IO_floatValue(line,positions,2_pInt) case ('divergence_correction') @@ -269,7 +282,7 @@ subroutine numerics_init case ('petsc_options') petsc_options = trim(line(positions(4):)) case ('myspectralsolver') - myspectralsolver = IO_stringValue(line,positions,2_pInt) + myspectralsolver = IO_lc(IO_stringValue(line,positions,2_pInt)) case ('err_f_tol') err_f_tol = IO_floatValue(line,positions,2_pInt) case ('err_p_tol')