diff --git a/code/DAMASK_spectral_Driver.f90 b/code/DAMASK_spectral_Driver.f90 index 720972d64..1ad0a5e25 100644 --- a/code/DAMASK_spectral_Driver.f90 +++ b/code/DAMASK_spectral_Driver.f90 @@ -37,9 +37,9 @@ program DAMASK_spectral_Driver use math use mesh, only : & - mesh_spectral_getResolution, & - mesh_spectral_getDimension, & - mesh_spectral_getHomogenization + res, & + geomdim, & + mesh_NcpElems use CPFEM, only: & CPFEM_initAll @@ -50,34 +50,34 @@ program DAMASK_spectral_Driver use numerics, only: & rotation_tol, & - myspectralsolver + mySpectralSolver use homogenization, only: & materialpoint_sizeResults, & materialpoint_results - - !use DAMASK_spectral_SolverAL + + use DAMASK_spectral_Utilities, only: & + boundaryCondition, & + solutionState, & + debugGeneral + use DAMASK_spectral_SolverBasic - use DAMASK_spectral_Utilities +!use DAMASK_spectral_SolverAL implicit none - type loadcase - real(pReal), dimension (3,3) :: deformation = 0.0_pReal, & ! applied velocity gradient or time derivative of deformation gradient - stress = 0.0_pReal, & ! stress BC (if applicable) - rotation = math_I3 ! rotation of BC (if applicable) + type loadCase + real(pReal), dimension (3,3) :: rotation = math_I3 ! rotation of BC + type(boundaryCondition) :: 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 integer(pInt) :: incs = 0_pInt, & ! number of increments outputfrequency = 1_pInt, & ! frequency of result writes restartfrequency = 0_pInt, & ! frequency of restart writes logscale = 0_pInt ! linear/logaritmic time inc flag - logical :: followFormerTrajectory = .true., & ! follow trajectory of former loadcase - velGradApplied = .false. ! decide wether velocity gradient or fdot is given - logical, dimension(3,3) :: maskDeformation = .false., & ! mask of deformation boundary conditions - maskStress = .false. ! mask of stress boundary conditions - logical, dimension(9) :: maskStressVector = .false. ! linear mask of boundary conditions - end type + logical :: followFormerTrajectory = .true. ! follow trajectory of former loadcase + end type loadCase !-------------------------------------------------------------------------------------------------- ! variables related to information from load case and geom file @@ -99,13 +99,11 @@ program DAMASK_spectral_Driver character(len=1024) :: & line - type(loadcase), allocatable, dimension(:) :: bc - type(solutionState) solres - type(BC_type) :: stress - + !-------------------------------------------------------------------------------------------------- ! loop variables, convergence etc. - 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 + 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 @@ -113,27 +111,31 @@ program DAMASK_spectral_Driver totalIncsCounter = 0_pInt,& notConvergedCounter = 0_pInt, convergedCounter = 0_pInt character(len=6) :: loadcase_string + + type(loadCase), allocatable, dimension(:) :: loadCases + type(solutionState) solres +!-------------------------------------------------------------------------------------------------- +! init DAMASK (all modules) call CPFEM_initAll(temperature = 300.0_pReal, element = 1_pInt, IP= 1_pInt) write(6,'(a)') '' write(6,'(a)') ' <<<+- DAMASK_spectral_Driver init -+>>>' write(6,'(a)') ' $Id$' #include "compilation_info.f90" - write(6,'(a,a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName()) - write(6,'(a,a)') ' Solver Job Name: ',trim(getSolverJobName()) + write(6,'(a,a)') ' Working Directory: ',trim(getSolverWorkingDirectoryName()) + write(6,'(a,a)') ' Solver Job Name: ',trim(getSolverJobName()) write(6,'(a)') '' write(6,'(a,a)') ' geometry file: ',trim(geometryFile) - write(6,'(a)') '=============================================================' - write(6,'(a,3(i12 ))') ' resolution a b c:', mesh_spectral_getResolution() - write(6,'(a,3(f12.5))') ' dimension x y z:', mesh_spectral_getDimension() - write(6,'(a,i5)') ' homogenization: ', mesh_spectral_getHomogenization() - write(6,'(a)') '=============================================================' - write(6,'(a,a)') 'Loadcase file: ',trim(loadCaseFile) + write(6,'(a)') '' + write(6,'(a,3(i12 ))') ' resolution a b c:', res + write(6,'(a,3(f12.5))') ' dimension x y z:', geomdim + write(6,'(a,i5)') ' homogenization: ', homog + write(6,'(a,a)') '','' + write(6,'(a,a)') ' Loadcase file: ',trim(loadCaseFile) write(6,'(a)') '' - !-------------------------------------------------------------------------------------------------- -! reading the load case file and allocate data structure containing load cases +! reading basic information from load case file and allocate data structure containing load cases call IO_open_file(myUnit,trim(loadCaseFile)) rewind(myUnit) do @@ -156,58 +158,62 @@ program DAMASK_spectral_Driver 100 if ((N_l + N_Fdot /= 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 (bc(N_n)) + allocate (loadCases(N_n)) + loadCases%P%myType='p' -!-------------------------------------------------------------------------------------------------- + !-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure rewind(myUnit) do read(myUnit,'(a1024)',END = 101) line if (IO_isBlank(line)) cycle ! skip empty lines - currentLoadcase = currentLoadcase + 1_pInt + currentLoadCase = currentLoadCase + 1_pInt positions = IO_stringPos(line,maxNchunksLoadcase) do i = 1_pInt,maxNchunksLoadcase select case (IO_lc(IO_stringValue(line,positions,i))) case('fdot','dotf','l','velocitygrad','velgrad','velocitygradient') ! assign values for the deformation BC matrix - bc(currentLoadcase)%velGradApplied = & - (IO_lc(IO_stringValue(line,positions,i)) == 'l'.or. & ! in case of given L, set flag to true - IO_lc(IO_stringValue(line,positions,i)) == 'velocitygrad'.or.& - IO_lc(IO_stringValue(line,positions,i)) == 'velgrad'.or.& - IO_lc(IO_stringValue(line,positions,i)) == 'velocitygradient') - temp_valueVector = 0.0_pReal - temp_maskVector = .false. + if (IO_lc(IO_stringValue(line,positions,i)) == 'l'.or. & ! in case of given L, set flag to true + IO_lc(IO_stringValue(line,positions,i)) == 'velocitygrad'.or.& + IO_lc(IO_stringValue(line,positions,i)) == 'velgrad'.or.& + IO_lc(IO_stringValue(line,positions,i)) == 'velocitygradient') then + loadCases(currentLoadCase)%deformation%myType = 'l' + else + loadCases(currentLoadCase)%deformation%myType = 'fdot' + endif forall (j = 1_pInt:9_pInt) temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' do j = 1_pInt,9_pInt if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j) enddo - bc(currentLoadcase)%maskDeformation = transpose(reshape(temp_maskVector,[ 3,3])) - bc(currentLoadcase)%deformation = math_plain9to33(temp_valueVector) + loadCases(currentLoadCase)%deformation%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) + loadCases(currentLoadCase)%deformation%maskFloat = merge(ones,zeroes,& + loadCases(currentLoadCase)%deformation%maskLogical) + loadCases(currentLoadCase)%deformation%values = math_plain9to33(temp_valueVector) case('p','pk1','piolakirchhoff','stress') temp_valueVector = 0.0_pReal - forall (j = 1_pInt:9_pInt) bc(currentLoadcase)%maskStressVector(j) =& - IO_stringValue(line,positions,i+j) /= '*' + forall (j = 1_pInt:9_pInt) temp_maskVector(j) = IO_stringValue(line,positions,i+j) /= '*' do j = 1_pInt,9_pInt - if (bc(currentLoadcase)%maskStressVector(j)) temp_valueVector(j) =& - IO_floatValue(line,positions,i+j) ! assign values for the bc(currentLoadcase)%stress matrix + if (temp_maskVector(j)) temp_valueVector(j) = IO_floatValue(line,positions,i+j) enddo - bc(currentLoadcase)%maskStress = transpose(reshape(bc(currentLoadcase)%maskStressVector,[ 3,3])) - bc(currentLoadcase)%stress = math_plain9to33(temp_valueVector) + loadCases(currentLoadCase)%P%maskLogical = transpose(reshape(temp_maskVector,[ 3,3])) + loadCases(currentLoadCase)%P%maskFloat = merge(ones,zeroes,& + loadCases(currentLoadCase)%P%maskLogical) + loadCases(currentLoadCase)%P%values = math_plain9to33(temp_valueVector) case('t','time','delta') ! increment time - bc(currentLoadcase)%time = IO_floatValue(line,positions,i+1_pInt) + loadCases(currentLoadCase)%time = IO_floatValue(line,positions,i+1_pInt) case('temp','temperature') ! starting temperature - bc(currentLoadcase)%temperature = IO_floatValue(line,positions,i+1_pInt) + loadCases(currentLoadCase)%temperature = IO_floatValue(line,positions,i+1_pInt) case('n','incs','increments','steps') ! number of increments - bc(currentLoadcase)%incs = IO_intValue(line,positions,i+1_pInt) + loadCases(currentLoadCase)%incs = IO_intValue(line,positions,i+1_pInt) case('logincs','logincrements','logsteps') ! number of increments (switch to log time scaling) - bc(currentLoadcase)%incs = IO_intValue(line,positions,i+1_pInt) - bc(currentLoadcase)%logscale = 1_pInt + loadCases(currentLoadCase)%incs = IO_intValue(line,positions,i+1_pInt) + loadCases(currentLoadCase)%logscale = 1_pInt case('f','freq','frequency','outputfreq') ! frequency of result writings - bc(currentLoadcase)%outputfrequency = IO_intValue(line,positions,i+1_pInt) + loadCases(currentLoadCase)%outputfrequency = IO_intValue(line,positions,i+1_pInt) case('r','restart','restartwrite') ! frequency of writing restart information - bc(currentLoadcase)%restartfrequency = max(0_pInt,IO_intValue(line,positions,i+1_pInt)) + loadCases(currentLoadCase)%restartfrequency = max(0_pInt,IO_intValue(line,positions,i+1_pInt)) case('guessreset','dropguessing') - bc(currentLoadcase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory - case('euler') ! rotation of currentLoadcase given in euler angles + loadCases(currentLoadCase)%followFormerTrajectory = .false. ! do not continue to predict deformation along former trajectory + case('euler') ! rotation of currentLoadCase given in euler angles l = 0_pInt ! assuming values given in radians k = 1_pInt ! assuming keyword indicating degree/radians select case (IO_lc(IO_stringValue(line,positions,i+1_pInt))) @@ -215,66 +221,75 @@ program DAMASK_spectral_Driver l = 1_pInt ! for conversion from degree to radian case('rad','radian') case default - k = 0_pInt ! immediately reading in angles, assuming radians + k = 0_pInt ! immediately readingk in angles, assuming radians end select forall(j = 1_pInt:3_pInt) temp33_Real(j,1) = & IO_floatValue(line,positions,i+k+j) * real(l,pReal) * inRad - bc(currentLoadcase)%rotation = math_EulerToR(temp33_Real(:,1)) - case('rotation','rot') ! assign values for the rotation of currentLoadcase matrix + loadCases(currentLoadCase)%rotation = math_EulerToR(temp33_Real(:,1)) + case('rotation','rot') ! assign values for the rotation of currentLoadCase matrix temp_valueVector = 0.0_pReal forall (j = 1_pInt:9_pInt) temp_valueVector(j) = IO_floatValue(line,positions,i+j) - bc(currentLoadcase)%rotation = math_plain9to33(temp_valueVector) + loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) end select enddo; enddo 101 close(myUnit) + !-------------------------------------------------------------------------------------------------- ! consistency checks and output of load case - bc(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadcase + loadCases(1)%followFormerTrajectory = .false. ! cannot guess along trajectory for first inc of first currentLoadCase errorID = 0_pInt - checkLoadcases: do currentLoadcase = 1_pInt, size(bc) - write (loadcase_string, '(i6)' ) currentLoadcase + checkLoadcases: do currentLoadCase = 1_pInt, size(loadCases) + write (loadcase_string, '(i6)' ) currentLoadCase - write(6,'(a)') '=============================================================' - write(6,'(a,i6)') 'currentLoadcase: ', currentLoadcase + write(6,'(2x,a,i6)') 'load case: ', currentLoadCase - if (.not. bc(currentLoadcase)%followFormerTrajectory) write(6,'(a)') 'drop guessing along trajectory' - if (bc(currentLoadcase)%velGradApplied) then + if (.not. loadCases(currentLoadCase)%followFormerTrajectory) write(6,'(2x,a)') 'drop guessing along trajectory' + if (loadCases(currentLoadCase)%deformation%myType=='l') then do j = 1_pInt, 3_pInt - if (any(bc(currentLoadcase)%maskDeformation(j,1:3) .eqv. .true.) .and. & - any(bc(currentLoadcase)%maskDeformation(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined + if (any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .true.) .and. & + any(loadCases(currentLoadCase)%deformation%maskLogical(j,1:3) .eqv. .false.)) errorID = 832_pInt ! each row should be either fully or not at all defined enddo - write(6,'(a)')'velocity gradient:' + write(6,'(2x,a)') 'velocity gradient:' else - write(6,'(a)')'deformation gradient rate:' + write(6,'(2x,a)') 'deformation gradient rate:' endif - write (6,'(3(3(f12.7,1x)/))',advance='no') merge(math_transpose33(bc(currentLoadcase)%deformation),& - reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(currentLoadcase)%maskDeformation)) - write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') ' stress / GPa:',& - 1e-9_pReal*merge(math_transpose33(bc(currentLoadcase)%stress),& - reshape(spread(DAMASK_NaN,1,9),[ 3,3]),transpose(bc(currentLoadcase)%maskStress)) - if (any(bc(currentLoadcase)%rotation /= math_I3)) & - write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') ' rotation of loadframe:',& - math_transpose33(bc(currentLoadcase)%rotation) - write(6,'(a,f12.6)') 'temperature:', bc(currentLoadcase)%temperature - write(6,'(a,f12.6)') 'time: ', bc(currentLoadcase)%time - write(6,'(a,i5)') 'increments: ', bc(currentLoadcase)%incs - write(6,'(a,i5)') 'output frequency: ', bc(currentLoadcase)%outputfrequency - write(6,'(a,i5)') 'restart frequency: ', bc(currentLoadcase)%restartfrequency + 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:',& + 1e-9_pReal*merge(math_transpose33(loadCases(currentLoadCase)%P%values),& + reshape(spread(DAMASK_NaN,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:',& + 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 - if (any(bc(currentLoadcase)%maskStress .eqv. bc(currentLoadcase)%maskDeformation)) errorID = 831_pInt ! exclusive or masking only - if (any(bc(currentLoadcase)%maskStress .and. transpose(bc(currentLoadcase)%maskStress) .and. & + 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]))) & errorID = 838_pInt ! no rotation is allowed by stress BC - if (any(abs(math_mul33x33(bc(currentLoadcase)%rotation,math_transpose33(bc(currentLoadcase)%rotation))& + if (any(abs(math_mul33x33(loadCases(currentLoadCase)%rotation,math_transpose33(loadCases(currentLoadCase)%rotation))& -math_I3) > reshape(spread(rotation_tol,1,9),[ 3,3]))& - .or. abs(math_det33(bc(currentLoadcase)%rotation)) > 1.0_pReal + rotation_tol)& + .or. abs(math_det33(loadCases(currentLoadCase)%rotation)) > 1.0_pReal + rotation_tol)& errorID = 846_pInt ! given rotation matrix contains strain - if (bc(currentLoadcase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment - if (bc(currentLoadcase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count - if (bc(currentLoadcase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency + if (loadCases(currentLoadCase)%time < 0.0_pReal) errorID = 834_pInt ! negative time increment + if (loadCases(currentLoadCase)%incs < 1_pInt) errorID = 835_pInt ! non-positive incs count + if (loadCases(currentLoadCase)%outputfrequency < 1_pInt) errorID = 836_pInt ! non-positive result frequency if (errorID > 0_pInt) call IO_error(error_ID = errorID, ext_msg = loadcase_string) enddo checkLoadcases + select case (myspectralsolver) + + case (DAMASK_spectral_SolverBasic_label) + call basic_init() + + !case (DAMASK_spectral_SolverAL_label) + ! call AL_init() + + end select !-------------------------------------------------------------------------------------------------- ! write header of output file if (appendToOutFile) then @@ -286,64 +301,53 @@ program DAMASK_spectral_Driver write(538) 'load', trim(loadCaseFile) write(538) 'workingdir', trim(getSolverWorkingDirectoryName()) write(538) 'geometry', trim(geometryFile) - write(538) 'resolution', mesh_spectral_getResolution() - write(538) 'dimension', mesh_spectral_getDimension() + write(538) 'resolution', res + write(538) 'dimension', geomdim write(538) 'materialpoint_sizeResults', materialpoint_sizeResults - write(538) 'loadcases', size(bc) - write(538) 'frequencies', bc%outputfrequency ! one entry per currentLoadcase - write(538) 'times', bc%time ! one entry per currentLoadcase - write(538) 'logscales', bc%logscale - write(538) 'increments', bc%incs ! one entry per currentLoadcase + write(538) 'loadcases', size(loadCases) + write(538) 'frequencies', loadCases%outputfrequency ! one entry per currentLoadCase + write(538) 'times', loadCases%time ! one entry per currentLoadCase + write(538) 'logscales', loadCases%logscale + write(538) 'increments', loadCases%incs ! one entry per currentLoadCase write(538) 'startingIncrement', restartInc - 1_pInt ! start with writing out the previous inc write(538) 'eoh' ! end of header - write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! initial (non-deformed or read-in) results + write(538) materialpoint_results(1_pInt:materialpoint_sizeResults,1,1_pInt:mesh_NcpElems) ! initial (non-deformed or read-in) results if (debugGeneral) write(6,'(a)') 'Header of result file written out' endif - select case (myspectralsolver) - - case (DAMASK_spectral_SolverBasic_label) - call basic_init() - - !case (DAMASK_spectral_SolverAL_label) - ! call AL_init() - - end select -!################################################################################################## -! Loop over loadcases defined in the currentLoadcase file -!################################################################################################## - loadCaseLooping: do currentLoadcase = 1_pInt, size(bc) - time0 = time ! currentLoadcase start time - if (bc(currentLoadcase)%followFormerTrajectory) then +!-------------------------------------------------------------------------------------------------- +! loopping over loadcases + loadCaseLooping: do currentLoadCase = 1_pInt, size(loadCases) + time0 = time ! currentLoadCase start time + if (loadCases(currentLoadCase)%followFormerTrajectory) then guessmode = 1.0_pReal else guessmode = 0.0_pReal ! change of load case, homogeneous guess for the first inc endif -!################################################################################################## -! loop oper incs defined in input file for current currentLoadcase -!################################################################################################## - incLooping: do inc = 1_pInt, bc(currentLoadcase)%incs +!-------------------------------------------------------------------------------------------------- +! loop oper incs defined in input file for current currentLoadCase + incLooping: do inc = 1_pInt, loadCases(currentLoadCase)%incs totalIncsCounter = totalIncsCounter + 1_pInt !-------------------------------------------------------------------------------------------------- ! forwarding time timeinc_old = timeinc - if (bc(currentLoadcase)%logscale == 0_pInt) then ! linear scale - timeinc = bc(currentLoadcase)%time/bc(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 (inc == 1_pInt) then ! 1st inc of 1st currentLoadcase of logarithmic scale - timeinc = bc(1)%time*(2.0_pReal**real( 1_pInt-bc(1)%incs ,pReal)) ! assume 1st inc is equal to 2nd - else ! not-1st inc of 1st currentLoadcase of logarithmic scale - timeinc = bc(1)%time*(2.0_pReal**real(inc-1_pInt-bc(1)%incs ,pReal)) + 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 + 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 - else ! not-1st currentLoadcase of logarithmic scale - timeinc = time0 *( (1.0_pReal + bc(currentLoadcase)%time/time0 )**(real( inc,pReal)/& - real(bc(currentLoadcase)%incs ,pReal))& - -(1.0_pReal + bc(currentLoadcase)%time/time0 )**(real( (inc-1_pInt),pReal)/& - real(bc(currentLoadcase)%incs ,pReal)) ) + else ! not-1st currentLoadCase of logarithmic scale + timeinc = time0 *( (1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( inc,pReal)/& + real(loadCases(currentLoadCase)%incs ,pReal))& + -(1.0_pReal + loadCases(currentLoadCase)%time/time0 )**(real( (inc-1_pInt),pReal)/& + real(loadCases(currentLoadCase)%incs ,pReal)) ) endif endif time = time + timeinc @@ -360,22 +364,20 @@ program DAMASK_spectral_Driver case (DAMASK_spectral_SolverBasic_label) solres = basic_solution (& guessmode,timeinc,timeinc_old, & - P_BC = bc(currentLoadcase)%stress, & - F_BC = bc(currentLoadcase)%deformation, & - ! temperature_bc = bc(currentLoadcase)%temperature, & - mask_stressVector = bc(currentLoadcase)%maskStressVector, & - velgrad = bc(currentLoadcase)%velGradApplied, & - rotation_BC = bc(currentLoadcase)%rotation) + 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 = bc(currentLoadcase)%stress, & - ! F_BC = bc(currentLoadcase)%deformation, & - ! ! temperature_bc = bc(currentLoadcase)%temperature, & - ! mask_stressVector = bc(currentLoadcase)%maskStressVector, & - ! velgrad = bc(currentLoadcase)%velGradApplied, & - ! rotation_BC = bc(currentLoadcase)%rotation) + ! P_BC = loadCases(currentLoadCase)%stress, & + ! F_BC = loadCases(currentLoadCase)%deformation, & + ! ! temperature_bc = loadCases(currentLoadCase)%temperature, & + ! mask_stressVector = loadCases(currentLoadCase)%maskStressVector, & + ! velgrad = loadCases(currentLoadCase)%velGradApplied, & + ! rotation_BC = loadCases(currentLoadCase)%rotation) end select @@ -389,17 +391,18 @@ program DAMASK_spectral_Driver notConvergedCounter = notConvergedCounter + 1_pInt endif - if (mod(inc,bc(currentLoadcase)%outputFrequency) == 0_pInt) then ! at output frequency + 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(1_pInt:materialpoint_sizeResults,1,1_pInt:Npoints) ! write result 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 + guessmode = 1.0_pReal ! keep guessing along former trajectory during same currentLoadCase enddo incLooping enddo loadCaseLooping + select case (myspectralsolver) case (DAMASK_spectral_SolverBasic_label) @@ -408,7 +411,8 @@ program DAMASK_spectral_Driver ! case (DAMASK_spectral_SolverAL_label) ! call AL_destroy() ! - end select + end select + write(6,'(a)') '' write(6,'(a)') '##################################################################' write(6,'(i6.6,a,i6.6,a,f5.1,a)') convergedCounter, ' out of ', & diff --git a/code/DAMASK_spectral_SolverAL.f90 b/code/DAMASK_spectral_SolverAL.f90 index 3f73ccc20..fdd2ce2f2 100644 --- a/code/DAMASK_spectral_SolverAL.f90 +++ b/code/DAMASK_spectral_SolverAL.f90 @@ -1,574 +1,590 @@ -module DAMASK_spectral_SolverAL +! module DAMASK_spectral_SolverAL - use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) + ! use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use DAMASK_spectral_Utilities + ! use DAMASK_spectral_Utilities - use math + ! use math - use mesh, only : & - mesh_spectral_getResolution, & - mesh_spectral_getDimension + ! use mesh, only : & + ! mesh_spectral_getResolution, & + ! mesh_spectral_getDimension - implicit none -#include -#include + ! implicit none +! #include +! #include - character (len=*), parameter, public :: & - DAMASK_spectral_SolverAL_label = 'AL' + ! character (len=*), parameter, public :: & + ! DAMASK_spectral_SolverAL_label = 'AL' -!-------------------------------------------------------------------------------------------------- -! PETSc data - SNES snes - KSP ksp - DM da - Vec x,r - PetscErrorCode ierr_psc - PetscMPIInt rank - PetscObject dummy - PetscInt xs,xm,gxs,gxm - PetscInt ys,ym,gys,gym - PetscInt zs,zm,gzs,gzm - character(len=1024) :: PetSc_options = '-snes_type ngmres -snes_ngmres_anderson -snes_monitor -snes_view' +! !-------------------------------------------------------------------------------------------------- +! ! PETSc data + ! SNES snes + ! KSP ksp + ! DM da + ! Vec x,r + ! PetscErrorCode ierr_psc + ! PetscMPIInt rank + ! PetscObject dummy + ! PetscInt xs,xm,gxs,gxm + ! PetscInt ys,ym,gys,gym + ! PetscInt zs,zm,gzs,gzm + ! character(len=1024) :: PetSc_options = '-snes_type ngmres -snes_ngmres_anderson -snes_monitor -snes_view' - external FormFunctionLocal, SNESConverged_Interactive + ! !external FormFunctionLocal, SNESConverged_Interactive -!-------------------------------------------------------------------------------------------------- -! common pointwise data - real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc, F_lambda, F_lambda_lastInc, P - real(pReal), dimension(:,:,:,:), allocatable :: coordinates - real(pReal), dimension(:,:,:), allocatable :: temperature +! !-------------------------------------------------------------------------------------------------- +! ! common pointwise data + ! real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc, F_lambda, F_lambda_lastInc, P + ! real(pReal), dimension(:,:,:,:), allocatable :: coordinates + ! real(pReal), dimension(:,:,:), allocatable :: temperature -!-------------------------------------------------------------------------------------------------- -! stress, stiffness and compliance average etc. - real(pReal), dimension(3,3) :: & - F_aim = math_I3, & - F_aim_lastInc = math_I3, & - P_av +! !-------------------------------------------------------------------------------------------------- +! ! stress, stiffness and compliance average etc. + ! real(pReal), dimension(3,3) :: & + ! F_aim = math_I3, & + ! F_aim_lastInc = math_I3, & + ! P_av - real(pReal), dimension(3,3,3,3) :: & - C_ref = 0.0_pReal, & - C = 0.0_pReal + ! real(pReal), dimension(3,3,3,3) :: & + ! !C_ref = 0.0_pReal, & + ! C = 0.0_pReal - integer(pInt) :: iter - real(pReal) :: err_div, err_stress + ! integer(pInt) :: iter + ! real(pReal) :: err_div, err_stress - contains + ! contains - subroutine AL_init() + ! subroutine AL_init() - 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 - implicit none - integer(pInt) :: i,j,k + ! implicit none + ! integer(pInt) :: i,j,k - call Utilities_init() + ! call Utilities_init() - allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lambda_lastInc(res(1),res(2),res(3),3,3), source = 0.0_pReal) - allocate (P ( res(1), res(2),res(3),3,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) + ! allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) + ! allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) + ! allocate (F_lambda ( res(1), res(2),res(3),3,3), source = 0.0_pReal) + ! allocate (F_lambda_lastInc(res(1),res(2),res(3),3,3), source = 0.0_pReal) + ! allocate (P ( res(1), res(2),res(3),3,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) -!-------------------------------------------------------------------------------------------------- -! init fields - if (restartInc == 1_pInt) then ! no deformation (no restart) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = math_I3 - F_lastInc(i,j,k,1:3,1:3) = math_I3 - F_lambda(i,j,k,1:3,1:3) = math_I3 - F_lambda_lastInc(i,j,k,1:3,1:3) = math_I3 - coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*res,pReal) - enddo; enddo; enddo - elseif (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& - restartInc - 1_pInt,' from file' - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& - trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& - trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',& - trim(getSolverJobName()),size(F_lambda)) - read (777,rec=1) F - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',& - trim(getSolverJobName()),size(F_lambda_lastInc)) - read (777,rec=1) F_lastInc - close (777) - call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_lastInc - close (777) +! !-------------------------------------------------------------------------------------------------- +! ! init fields + ! if (restartInc == 1_pInt) then ! no deformation (no restart) + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! F(i,j,k,1:3,1:3) = math_I3 + ! F_lastInc(i,j,k,1:3,1:3) = math_I3 + ! F_lambda(i,j,k,1:3,1:3) = math_I3 + ! F_lambda_lastInc(i,j,k,1:3,1:3) = math_I3 + ! coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & + ! - geomdim/real(2_pInt*res,pReal) + ! enddo; enddo; enddo + ! elseif (restartInc > 1_pInt) then ! using old values from file + ! if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& + ! restartInc - 1_pInt,' from file' + ! call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& + ! trim(getSolverJobName()),size(F)) + ! read (777,rec=1) F + ! close (777) + ! call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& + ! trim(getSolverJobName()),size(F_lastInc)) + ! read (777,rec=1) F_lastInc + ! close (777) + ! call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda',& + ! trim(getSolverJobName()),size(F_lambda)) + ! read (777,rec=1) F + ! close (777) + ! call IO_read_jobBinaryFile(777,'convergedSpectralDefgradLambda_lastInc',& + ! trim(getSolverJobName()),size(F_lambda_lastInc)) + ! read (777,rec=1) F_lastInc + ! close (777) + ! call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) + ! read (777,rec=1) F_aim + ! close (777) + ! call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) + ! read (777,rec=1) F_aim_lastInc + ! close (777) - coordinates = 0.0 ! change it later!!! - endif + ! coordinates = 0.0 ! change it later!!! + ! endif - call constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& - P,C,P_av,.false.,math_I3) + ! call constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& + ! P,C,P_av,.false.,math_I3) -!-------------------------------------------------------------------------------------------------- -! reference stiffness - if (restartInc == 1_pInt) then - call IO_write_jobBinaryFile(777,'C_ref',size(C)) - write (777,rec=1) C - close(777) - elseif (restartInc > 1_pInt) then - call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) - read (777,rec=1) C - close (777) - endif +! !-------------------------------------------------------------------------------------------------- +! ! reference stiffness + ! if (restartInc == 1_pInt) then + ! call IO_write_jobBinaryFile(777,'C_ref',size(C)) + ! write (777,rec=1) C + ! close(777) + ! elseif (restartInc > 1_pInt) then + ! call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) + ! read (777,rec=1) C + ! close (777) + ! endif - call Utilities_updateGamma(C_ref) + ! call Utilities_updateGamma(C_ref) -!-------------------------------------------------------------------------------------------------- -! PETSc Init - call PetscInitialize(PETSC_NULL_CHARACTER,ierr_psc) - call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr_psc) +! !-------------------------------------------------------------------------------------------------- +! ! PETSc Init + ! call PetscInitialize(PETSC_NULL_CHARACTER,ierr_psc) + ! call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr_psc) - call SNESCreate(PETSC_COMM_WORLD,snes,ierr_psc) - call DMDACreate3d(PETSC_COMM_WORLD, & - DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & - DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & - 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr_psc) - call DMCreateGlobalVector(da,x,ierr_psc) - call VecDuplicate(x,r,ierr_psc) - call DMDASetLocalFunction(da,FormFunctionLocal,ierr_psc) + ! call SNESCreate(PETSC_COMM_WORLD,snes,ierr_psc) + ! call DMDACreate3d(PETSC_COMM_WORLD, & + ! DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE, & + ! DMDA_STENCIL_BOX,res(1),res(2),res(3),PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE, & + ! 18,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr_psc) + ! call DMCreateGlobalVector(da,x,ierr_psc) + ! call VecDuplicate(x,r,ierr_psc) + ! call DMDASetLocalFunction(da,FormFunctionLocal,ierr_psc) - call SNESSetDM(snes,da,ierr_psc) - call SNESSetFunction(snes,r,SNESDMDAComputeFunction,da,ierr_psc) - call SNESSetConvergenceTest(snes,SNESConverged_Interactive,dummy,PETSC_NULL_FUNCTION,ierr_psc) - call PetscOptionsInsertString(PetSc_options,ierr_psc) - call SNESSetFromOptions(snes,ierr_psc) - call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr_psc) - call DMDAGetCorners(da,gxs,gys,gzs,gxm,gym,gzm,ierr_psc) + ! call SNESSetDM(snes,da,ierr_psc) + ! call SNESSetFunction(snes,r,SNESDMDAComputeFunction,da,ierr_psc) + ! call SNESSetConvergenceTest(snes,SNESConverged_Interactive,dummy,PETSC_NULL_FUNCTION,ierr_psc) + ! call PetscOptionsInsertString(PetSc_options,ierr_psc) + ! call SNESSetFromOptions(snes,ierr_psc) + ! call DMDAGetCorners(da,xs,ys,zs,xm,ym,zm,ierr_psc) + ! call DMDAGetCorners(da,gxs,gys,gzs,gxm,gym,gzm,ierr_psc) - xs = xs+1; gxs = gxs+1; xm = xm-1; gxm = gxm-1 - ys = ys+1; gys = gys+1; ym = ym-1; gym = gym-1 - zs = zs+1; gzs = gzs+1; zm = zm-1; gzm = gzm-1 + ! xs = xs+1; gxs = gxs+1; xm = xm-1; gxm = gxm-1 + ! ys = ys+1; gys = gys+1; ym = ym-1; gym = gym-1 + ! zs = zs+1; gzs = gzs+1; zm = zm-1; gzm = gzm-1 - end subroutine AL_init + ! end subroutine AL_init -type(solutionState) function AL_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,mask_stressVector,velgrad,rotation_BC) +! type(solutionState) function AL_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,mask_stressVector,velgrad,rotation_BC) - use numerics, only: & - itmax, & - itmin, & - update_gamma + ! use numerics, only: & + ! itmax, & + ! itmin, & + ! update_gamma - use IO, only: & - IO_write_JobBinaryFile + ! use IO, only: & + ! IO_write_JobBinaryFile - use FEsolving, only: & - restartWrite + ! use FEsolving, only: & + ! restartWrite - implicit none + ! implicit none -!-------------------------------------------------------------------------------------------------- -! input data for solution +! !-------------------------------------------------------------------------------------------------- +! ! input data for solution - real(pReal), intent(in) :: timeinc, timeinc_old - real(pReal), intent(in) :: guessmode - logical, intent(in) :: velgrad - real(pReal), dimension(3,3), intent(in) :: P_BC,F_BC,rotation_BC - logical, dimension(9), intent(in) :: mask_stressVector + ! real(pReal), intent(in) :: timeinc, timeinc_old + ! real(pReal), intent(in) :: guessmode + ! logical, intent(in) :: velgrad + ! real(pReal), dimension(3,3), intent(in) :: P_BC,F_BC,rotation_BC + ! logical, dimension(9), intent(in) :: mask_stressVector -!-------------------------------------------------------------------------------------------------- -! loop variables, convergence etc. +! !-------------------------------------------------------------------------------------------------- +! ! loop variables, convergence etc. - real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal - real(pReal), dimension(3,3) :: temp33_Real - real(pReal), dimension(3,3,3,3) :: S - real(pReal), dimension(3,3) :: mask_stress, & - mask_defgrad, & - deltaF_aim, & - F_aim_lab, & - F_aim_lab_lastIter - integer(pInt) :: i, j, k - logical :: ForwardData - real(pReal) :: defgradDet - real(pReal) :: defgradDetMax, defgradDetMin + ! real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal + ! real(pReal), dimension(3,3) :: temp33_Real + ! real(pReal), dimension(3,3,3,3) :: S + ! real(pReal), dimension(3,3) :: mask_stress, & + ! mask_defgrad, & + ! deltaF_aim, & + ! F_aim_lab, & + ! F_aim_lab_lastIter + ! integer(pInt) :: i, j, k + ! logical :: ForwardData + ! real(pReal) :: defgradDet + ! real(pReal) :: defgradDetMax, defgradDetMin - PetscScalar, pointer :: xx_psc(:) + ! PetscScalar, pointer :: xx_psc(:) - mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) - mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) + ! mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) + ! mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) - if (restartWrite) then - write(6,'(a)') 'writing converged results for restart' - call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) ! writing deformation gradient field to file - write (777,rec=1) F_LastInc - close (777) - call IO_write_jobBinaryFile(777,'C',size(C)) - write (777,rec=1) C - close(777) - endif + ! if (restartWrite) then + ! write(6,'(a)') 'writing converged results for restart' + ! call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) ! writing deformation gradient field to file + ! write (777,rec=1) F_LastInc + ! close (777) + ! call IO_write_jobBinaryFile(777,'C',size(C)) + ! write (777,rec=1) C + ! close(777) + ! endif - ForwardData = .True. - if (velgrad) then ! calculate deltaF_aim from given L and current F - deltaF_aim = timeinc * mask_defgrad * math_mul33x33(F_BC, F_aim) - else ! deltaF_aim = fDot *timeinc where applicable - deltaF_aim = timeinc * mask_defgrad * F_BC - endif + ! ForwardData = .True. + ! if (velgrad) then ! calculate deltaF_aim from given L and current F + ! deltaF_aim = timeinc * mask_defgrad * math_mul33x33(F_BC, F_aim) + ! else ! deltaF_aim = fDot *timeinc where applicable + ! deltaF_aim = timeinc * mask_defgrad * F_BC + ! endif -!-------------------------------------------------------------------------------------------------- -! winding forward of deformation aim in loadcase system - temp33_Real = F_aim - F_aim = F_aim & - + guessmode * mask_stress * (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 +! !-------------------------------------------------------------------------------------------------- +! ! winding forward of deformation aim in loadcase system + ! temp33_Real = F_aim + ! F_aim = F_aim & + ! + guessmode * mask_stress * (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) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_Real = F(i,j,k,1:3,1:3) - F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon - + guessmode * (F(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))* & - timeinc/timeinc_old + (1.0_pReal-guessmode) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable - F_lastInc(i,j,k,1:3,1:3) = temp33_Real - temp33_Real = F_lambda(i,j,k,1:3,1:3) - F_lambda(i,j,k,1:3,1:3) = F_lambda(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon - + guessmode * (F_lambda(i,j,k,1:3,1:3) - F_lambda_lastInc(i,j,k,1:3,1:3))* & - timeinc/timeinc_old + (1.0_pReal-guessmode) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable - F_lambda_lastInc(i,j,k,1:3,1:3) = temp33_Real - enddo; enddo; enddo - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),& ! calculate current coordinates - 1.0_pReal,F_lastInc,coordinates) +! !-------------------------------------------------------------------------------------------------- +! ! update local deformation gradient and coordinates + ! deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! temp33_Real = F(i,j,k,1:3,1:3) + ! F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon + ! + guessmode * (F(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))* & + ! timeinc/timeinc_old + (1.0_pReal-guessmode) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable + ! F_lastInc(i,j,k,1:3,1:3) = temp33_Real + ! temp33_Real = F_lambda(i,j,k,1:3,1:3) + ! F_lambda(i,j,k,1:3,1:3) = F_lambda(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon + ! + guessmode * (F_lambda(i,j,k,1:3,1:3) - F_lambda_lastInc(i,j,k,1:3,1:3))* & + ! timeinc/timeinc_old + (1.0_pReal-guessmode) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable + ! F_lambda_lastInc(i,j,k,1:3,1:3) = temp33_Real + ! enddo; enddo; enddo + ! call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),& ! calculate current coordinates + ! 1.0_pReal,F_lastInc,coordinates) - iter = 0_pInt - S = Utilities_stressBC(rotation_BC,mask_stressVector,C) - if (update_gamma) call Utilities_updateGamma(C) + ! iter = 0_pInt + ! S = Utilities_stressBC(rotation_BC,mask_stressVector,C) + ! if (update_gamma) call Utilities_updateGamma(C) - call VecGetArrayF90(x,xx_psc,ierr_psc) - call FormInitialGuessLocal(xx_psc) - call VecRestoreArrayF90(x,xx_psc,ierr_psc) - call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr_psc) + ! call VecGetArrayF90(x,xx_psc,ierr_psc) + ! call FormInitialGuessLocal(xx_psc) + ! call VecRestoreArrayF90(x,xx_psc,ierr_psc) + ! call SNESSolve(snes,PETSC_NULL_OBJECT,x,ierr_psc) - convergenceLoop: do while((iter < itmax .and. (any([err_div ,err_stress] > 1.0_pReal)))& - .or. iter < itmin) + ! convergenceLoop: do while((iter < itmax .and. (any([err_div ,err_stress] > 1.0_pReal)))& + ! .or. iter < itmin) - iter = iter + 1_pInt -!-------------------------------------------------------------------------------------------------- -! report begin of new iteration - write(6,'(a)') '' - write(6,'(a)') '==================================================================' - write(6,'(3(a,i6.6))') ' @ 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) + ! iter = iter + 1_pInt +! !-------------------------------------------------------------------------------------------------- +! ! report begin of new iteration + ! write(6,'(a)') '' + ! write(6,'(a)') '==================================================================' + ! write(6,'(3(a,i6.6))') ' @ 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 - call constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,& - P,C,P_av,ForwardData,rotation_BC) - ForwardData = .False. +! !-------------------------------------------------------------------------------------------------- +! ! evaluate constitutive response + ! call constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,& + ! P,C,P_av,ForwardData,rotation_BC) + ! ForwardData = .False. -!-------------------------------------------------------------------------------------------------- -! stress BC handling - if(any(mask_stressVector)) then ! calculate stress BC if applied - F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC))) - err_stress = mask_stress * (P_av - P_BC))) - else - err_stress = 0.0_pReal - endif +! !-------------------------------------------------------------------------------------------------- +! ! stress BC handling + ! if(any(mask_stressVector)) then ! calculate stress BC if applied + ! F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC))) + ! err_stress = mask_stress * (P_av - P_BC) + ! else + ! err_stress = 0.0_pReal + ! endif - F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame + ! F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) ! boundary conditions from load frame into lab (Fourier) frame -!-------------------------------------------------------------------------------------------------- -! updated deformation gradient - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = P - call FFT_forward() - err_div = calcDivergence() - call convolution_fourier(F_aim_lab_lastIter - F_aim_lab, C_ref) - call FFT_backward() +! !-------------------------------------------------------------------------------------------------- +! ! updated deformation gradient + ! field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = P + ! call FFT_forward() + ! err_div = calcDivergence() + ! call convolution_fourier(F_aim_lab_lastIter - F_aim_lab, C_ref) + ! call FFT_backward() - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization - enddo; enddo; enddo + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization + ! enddo; enddo; enddo -!-------------------------------------------------------------------------------------------------- -! calculate some additional output - if(debugGeneral) then - maxCorrectionSkew = 0.0_pReal - maxCorrectionSym = 0.0_pReal - temp33_Real = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - maxCorrectionSym = max(maxCorrectionSym,& - maxval(math_symmetric33(field_real(i,j,k,1:3,1:3)))) - maxCorrectionSkew = max(maxCorrectionSkew,& - maxval(math_skew33(field_real(i,j,k,1:3,1:3)))) - temp33_Real = temp33_Real + field_real(i,j,k,1:3,1:3) - enddo; enddo; enddo - write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',& - maxCorrectionSym*wgt - write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',& - maxCorrectionSkew*wgt - write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',& - maxval(math_symmetric33(temp33_real))/& - maxval(math_skew33(temp33_real)) - endif +! !-------------------------------------------------------------------------------------------------- +! ! calculate some additional output + ! if(debugGeneral) then + ! maxCorrectionSkew = 0.0_pReal + ! maxCorrectionSym = 0.0_pReal + ! temp33_Real = 0.0_pReal + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! maxCorrectionSym = max(maxCorrectionSym,& + ! maxval(math_symmetric33(field_real(i,j,k,1:3,1:3)))) + ! maxCorrectionSkew = max(maxCorrectionSkew,& + ! maxval(math_skew33(field_real(i,j,k,1:3,1:3)))) + ! temp33_Real = temp33_Real + field_real(i,j,k,1:3,1:3) + ! enddo; enddo; enddo + ! write(6,'(a,1x,es11.4)') 'max symmetric correction of deformation =',& + ! maxCorrectionSym*wgt + ! write(6,'(a,1x,es11.4)') 'max skew correction of deformation =',& + ! maxCorrectionSkew*wgt + ! write(6,'(a,1x,es11.4)') 'max sym/skew of avg correction = ',& + ! maxval(math_symmetric33(temp33_real))/& + ! maxval(math_skew33(temp33_real)) + ! endif -!-------------------------------------------------------------------------------------------------- -! calculate bounds of det(F) and report - if(debugGeneral) then - defgradDetMax = -huge(1.0_pReal) - defgradDetMin = +huge(1.0_pReal) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - defgradDet = math_det33(F(i,j,k,1:3,1:3)) - defgradDetMax = max(defgradDetMax,defgradDet) - defgradDetMin = min(defgradDetMin,defgradDet) - enddo; enddo; enddo +! !-------------------------------------------------------------------------------------------------- +! ! calculate bounds of det(F) and report + ! if(debugGeneral) then + ! defgradDetMax = -huge(1.0_pReal) + ! defgradDetMin = +huge(1.0_pReal) + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! defgradDet = math_det33(F(i,j,k,1:3,1:3)) + ! defgradDetMax = max(defgradDetMax,defgradDet) + ! defgradDetMin = min(defgradDetMin,defgradDet) + ! enddo; enddo; enddo - write(6,'(a,1x,es11.4)') 'max determinant of deformation =', defgradDetMax - write(6,'(a,1x,es11.4)') 'min determinant of deformation =', defgradDetMin - endif - enddo convergenceLoop + ! write(6,'(a,1x,es11.4)') 'max determinant of deformation =', defgradDetMax + ! write(6,'(a,1x,es11.4)') 'min determinant of deformation =', defgradDetMin + ! endif + ! enddo convergenceLoop -end function AL_solution +! end function AL_solution -subroutine AL_destroy() +! subroutine AL_destroy() -implicit none +! implicit none -call VecDestroy(x,ierr_psc) -call VecDestroy(r,ierr_psc) -call SNESDestroy(snes,ierr_psc) -call DMDestroy(da,ierr_psc) -call PetscFinalize(ierr_psc) -call Utilities_destroy() +! call VecDestroy(x,ierr_psc) +! call VecDestroy(r,ierr_psc) +! call SNESDestroy(snes,ierr_psc) +! call DMDestroy(da,ierr_psc) +! call PetscFinalize(ierr_psc) +! call Utilities_destroy() -end subroutine AL_destroy +! end subroutine AL_destroy -! ------------------------------------------------------------------- +! ! ------------------------------------------------------------------- -subroutine FormInitialGuessLocal(xx_psc) +! subroutine FormInitialGuessLocal(xx_psc) - implicit none -#include + ! implicit none +! #include -! Input/output variables: +! ! Input/output variables: - PetscScalar xx_psc(0:17,gxs:(gxs+gxm),gys:(gys+gym),gxs:(gzs+gzm)) - integer(pInt) :: i, j, k + ! PetscScalar xx_psc(0:17,gxs:(gxs+gxm),gys:(gys+gym),gxs:(gzs+gzm)) + ! integer(pInt) :: i, j, k -! Compute function over the locally owned part of the grid +! ! Compute function over the locally owned part of the grid - do k=gzs,gzs+gzm; do j=gys,gys+gym; do i=gxs,gxs+gxm - xx_psc(0,i,j,k) = F(i,j,k,1,1) - xx_psc(1,i,j,k) = F(i,j,k,1,2) - xx_psc(2,i,j,k) = F(i,j,k,1,3) - xx_psc(3,i,j,k) = F(i,j,k,2,1) - xx_psc(4,i,j,k) = F(i,j,k,2,2) - xx_psc(5,i,j,k) = F(i,j,k,2,3) - xx_psc(6,i,j,k) = F(i,j,k,3,1) - xx_psc(7,i,j,k) = F(i,j,k,3,2) - xx_psc(8,i,j,k) = F(i,j,k,3,3) - xx_psc(9,i,j,k) = F_lambda(i,j,k,1,1) - xx_psc(10,i,j,k) = F_lambda(i,j,k,1,2) - xx_psc(11,i,j,k) = F_lambda(i,j,k,1,3) - xx_psc(12,i,j,k) = F_lambda(i,j,k,2,1) - xx_psc(13,i,j,k) = F_lambda(i,j,k,2,2) - xx_psc(14,i,j,k) = F_lambda(i,j,k,2,3) - xx_psc(15,i,j,k) = F_lambda(i,j,k,3,1) - xx_psc(16,i,j,k) = F_lambda(i,j,k,3,2) - xx_psc(17,i,j,k) = F_lambda(i,j,k,3,3) - enddo; enddo; enddo + ! do k=gzs,gzs+gzm; do j=gys,gys+gym; do i=gxs,gxs+gxm + ! xx_psc(0,i,j,k) = F(i,j,k,1,1) + ! xx_psc(1,i,j,k) = F(i,j,k,1,2) + ! xx_psc(2,i,j,k) = F(i,j,k,1,3) + ! xx_psc(3,i,j,k) = F(i,j,k,2,1) + ! xx_psc(4,i,j,k) = F(i,j,k,2,2) + ! xx_psc(5,i,j,k) = F(i,j,k,2,3) + ! xx_psc(6,i,j,k) = F(i,j,k,3,1) + ! xx_psc(7,i,j,k) = F(i,j,k,3,2) + ! xx_psc(8,i,j,k) = F(i,j,k,3,3) + ! xx_psc(9,i,j,k) = F_lambda(i,j,k,1,1) + ! xx_psc(10,i,j,k) = F_lambda(i,j,k,1,2) + ! xx_psc(11,i,j,k) = F_lambda(i,j,k,1,3) + ! xx_psc(12,i,j,k) = F_lambda(i,j,k,2,1) + ! xx_psc(13,i,j,k) = F_lambda(i,j,k,2,2) + ! xx_psc(14,i,j,k) = F_lambda(i,j,k,2,3) + ! xx_psc(15,i,j,k) = F_lambda(i,j,k,3,1) + ! xx_psc(16,i,j,k) = F_lambda(i,j,k,3,2) + ! xx_psc(17,i,j,k) = F_lambda(i,j,k,3,3) + ! enddo; enddo; enddo - return -end subroutine FormInitialGuessLocal + ! return +! end subroutine FormInitialGuessLocal -! --------------------------------------------------------------------- -! -! Input Parameter: -! x - local vector data -! -! Output Parameters: -! f - local vector data, f(x) -! ierr - error code -! -! Notes: -! This routine uses standard Fortran-style computations over a 3-dim array. -! -subroutine FormFunctionLocal(in,x_scal,f_scal,dummy,ierr_psc) +! ! --------------------------------------------------------------------- +! ! +! ! Input Parameter: +! ! x - local vector data +! ! +! ! Output Parameters: +! ! f - local vector data, f(x) +! ! ierr - error code +! ! +! ! Notes: +! ! This routine uses standard Fortran-style computations over a 3-dim array. +! ! +! subroutine FormFunctionLocal(in,x_scal,f_scal,dummy,ierr_psc) - use numerics, only: & - itmax, & - itmin + ! use numerics, only: & + ! itmax, & + ! itmin - implicit none -#include + ! implicit none +! #include -! Input/output variables: - DMDALocalInfo in(DMDA_LOCAL_INFO_SIZE) - PetscScalar x_scal(0:17,XG_RANGE,YG_RANGE,ZG_RANGE) - PetscScalar f_scal(0:17,X_RANGE,Y_RANGE,Z_RANGE) - real(pReal), dimension (3,3) :: temp - PetscObject dummy - -! Compute function over the locally owned part of the grid +! integer(pInt) :: i,j,k +! ! Input/output variables: + ! DMDALocalInfo in(DMDA_LOCAL_INFO_SIZE) + ! PetscScalar x_scal(0:17,XG_RANGE,YG_RANGE,ZG_RANGE) + ! PetscScalar f_scal(0:17,X_RANGE,Y_RANGE,Z_RANGE) + ! real(pReal), dimension (3,3) :: temp, lambda_av, F_star_av + ! PetscObject dummy + ! PetscInt gzs,gze + ! PetscInt gxs,gxe + ! PetscInt gys,gye + ! PetscErrorCode ierr_psc - iter = iter + 1_pInt + + ! gxs = in(DMDA_LOCAL_INFO_GXS)+1 + ! gxe = gxs+in(DMDA_LOCAL_INFO_GXM)-1 + ! gys = in(DMDA_LOCAL_INFO_GYS)+1 + ! gye = gxs+in(DMDA_LOCAL_INFO_GYM)-1 + ! gzs = in(DMDA_LOCAL_INFO_GZS)+1 + ! gze = gzs+in(DMDA_LOCAL_INFO_GZM)-1 + + ! iter = iter + 1_pInt -!-------------------------------------------------------------------------------------------------- -! report begin of new iteration - write(6,'(a)') '' - write(6,'(a)') '==================================================================' - write(6,'(3(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax - write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& - math_transpose33(F_aim) +! !-------------------------------------------------------------------------------------------------- +! ! report begin of new iteration + ! write(6,'(a)') '' + ! write(6,'(a)') '==================================================================' + ! write(6,'(3(a,i6.6))') ' @ Iter. ',itmin,' < ',iter,' < ',itmax + ! write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& + ! math_transpose33(F_aim) - F_star_av = 0.0 - lambda_av = 0.0 - do k=gzs,gze; do j=gys,gye; do i=gxs,gxe - F(i,j,k,1,1) = x_scal(0,i,j,k) - F(i,j,k,1,2) = x_scal(1,i,j,k) - F(i,j,k,1,3) = x_scal(2,i,j,k) - F(i,j,k,2,1) = x_scal(3,i,j,k) - F(i,j,k,2,2) = x_scal(4,i,j,k) - F(i,j,k,2,3) = x_scal(5,i,j,k) - F(i,j,k,3,1) = x_scal(6,i,j,k) - F(i,j,k,3,2) = x_scal(7,i,j,k) - F(i,j,k,3,3) = x_scal(8,i,j,k) - F_lambda(i,j,k,1,1) = x_scal(9,i,j,k) - F_lambda(i,j,k,1,2) = x_scal(10,i,j,k) - F_lambda(i,j,k,1,3) = x_scal(11,i,j,k) - F_lambda(i,j,k,2,1) = x_scal(12,i,j,k) - F_lambda(i,j,k,2,2) = x_scal(13,i,j,k) - F_lambda(i,j,k,2,3) = x_scal(14,i,j,k) - F_lambda(i,j,k,3,1) = x_scal(15,i,j,k) - F_lambda(i,j,k,3,2) = x_scal(16,i,j,k) - F_lambda(i,j,k,3,3) = x_scal(17,i,j,k) - F_star_av = F_star_av + F(i,j,k,1:3,1:3) - lambda_av = lambda_av + F_lambda(i,j,k,1:3,1:3) - enddo; enddo; enddo - F_star_av = F_star_av *wgt - lambda_av = math_mul3333xx33(C_inc0,lambda_av*wgt-math_I3) +! ! F_star_av = 0.0 +! ! lambda_av = 0.0 +! ! do k=gzs,gze; do j=gys,gye; do i=gxs,gxe +! ! F(i,j,k,1,1) = x_scal(0,i,j,k) +! ! F(i,j,k,1,2) = x_scal(1,i,j,k) +! ! F(i,j,k,1,3) = x_scal(2,i,j,k) +! ! F(i,j,k,2,1) = x_scal(3,i,j,k) +! ! F(i,j,k,2,2) = x_scal(4,i,j,k) +! ! F(i,j,k,2,3) = x_scal(5,i,j,k) +! ! F(i,j,k,3,1) = x_scal(6,i,j,k) +! ! F(i,j,k,3,2) = x_scal(7,i,j,k) +! ! F(i,j,k,3,3) = x_scal(8,i,j,k) +! ! F_lambda(i,j,k,1,1) = x_scal(9,i,j,k) +! ! F_lambda(i,j,k,1,2) = x_scal(10,i,j,k) +! ! F_lambda(i,j,k,1,3) = x_scal(11,i,j,k) +! ! F_lambda(i,j,k,2,1) = x_scal(12,i,j,k) +! ! F_lambda(i,j,k,2,2) = x_scal(13,i,j,k) +! ! F_lambda(i,j,k,2,3) = x_scal(14,i,j,k) +! ! F_lambda(i,j,k,3,1) = x_scal(15,i,j,k) +! ! F_lambda(i,j,k,3,2) = x_scal(16,i,j,k) +! ! F_lambda(i,j,k,3,3) = x_scal(17,i,j,k) +! ! F_star_av = F_star_av + F(i,j,k,1:3,1:3) +! ! lambda_av = lambda_av + F_lambda(i,j,k,1:3,1:3) +! ! enddo; enddo; enddo +! ! F_star_av = F_star_av *wgt + ! lambda_av = math_mul3333xx33(C_inc0,lambda_av*wgt-math_I3) -!-------------------------------------------------------------------------------------------------- -! evaluate constitutive response - call constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,& - P,C,P_av,ForwardData,rotation_BC) - ForwardData = .False. +! !-------------------------------------------------------------------------------------------------- +! ! evaluate constitutive response + ! real(pReal), intent(in) :: timeinc, timeinc_old + ! real(pReal), intent(in) :: guessmode + ! logical, intent(in) :: velgrad + ! real(pReal), dimension(3,3), intent(in) :: P_BC,F_BC,rotation_BC + ! logical, dimension(9), intent(in) :: mask_stressVector + ! call constitutiveResponse(coordinates,F,F_lastInc,temperature,timeinc,& + ! P,C,P_av,ForwardData,rotation_BC) + ! ForwardData = .False. -!-------------------------------------------------------------------------------------------------- -! stress BC handling - if(any(mask_stressVector)) then ! calculate stress BC if applied - F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC))) - err_stress = mask_stress * (P_av - P_BC))) - else - err_stress = 0.0_pReal - endif +! !-------------------------------------------------------------------------------------------------- +! ! stress BC handling + ! if(any(mask_stressVector)) then ! calculate stress BC if applied + ! F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC))) + ! err_stress = maxval(mask_stress * (P_av - P_BC)) + ! else + ! err_stress = 0.0_pReal + ! endif - F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) + ! F_aim_lab = math_rotate_backward33(F_aim,rotation_BC) -!-------------------------------------------------------------------------------------------------- -! doing Fourier transform - field_real = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_ref,F_lambda(i,j,k,1:3,1:3)-F(i,j,k,1:3,1:3)) +! !-------------------------------------------------------------------------------------------------- +! ! doing Fourier transform + ! field_real = 0.0_pReal + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! field_real(i,j,k,1:3,1:3) = math_mul3333xx33(C_ref,F_lambda(i,j,k,1:3,1:3)-F(i,j,k,1:3,1:3)) - enddo; enddo; enddo + ! enddo; enddo; enddo - call Utilities_forwardFFT() - call Utilities_fourierConvolution(F_aim_lab) - call Utilities_backwardFFT() + ! call Utilities_forwardFFT() + ! call Utilities_fourierConvolution(F_aim_lab) + ! call Utilities_backwardFFT() - err_f = 0.0_pReal - err_f_point = 0.0_pReal - err_p = 0.0_pReal - err_p_point = 0.0_pReal + ! err_f = 0.0_pReal + ! err_f_point = 0.0_pReal + ! err_p = 0.0_pReal + ! err_p_point = 0.0_pReal - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_real = field_real(i,j,k,1:3,1:3) - F(i,j,k,1:3,1:3) - err_f_point = max(err_f_point, maxval(abs(temp33_real))) - err_f = err_f + sum(temp33_real*temp33_real) + ! do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + ! temp33_real = field_real(i,j,k,1:3,1:3) - F(i,j,k,1:3,1:3) + ! err_f_point = max(err_f_point, maxval(abs(temp33_real))) + ! err_f = err_f + sum(temp33_real*temp33_real) - temp33_real = F_lambda(i,j,k,1:3,1:3) - & - math_mul3333xx33(S_inc0,P(i,j,k,1:3,1:3)) + math_I3 - err_p_point = max(err_p_point, maxval(abs(temp33_real))) - err_p = err_p + sum(temp33_real*temp33_real) - enddo; enddo; enddo + ! temp33_real = F_lambda(i,j,k,1:3,1:3) - & + ! math_mul3333xx33(S_inc0,P(i,j,k,1:3,1:3)) + math_I3 + ! err_p_point = max(err_p_point, maxval(abs(temp33_real))) + ! err_p = err_p + sum(temp33_real*temp33_real) + ! enddo; enddo; enddo - err_f = wgt*sqrt(err_f/sum((F_aim-math_I3)*(F_aim-math_I3))) - err_p = wgt*sqrt(err_p/sum((F_aim-math_I3)*(F_aim-math_I3))) + ! err_f = wgt*sqrt(err_f/sum((F_aim-math_I3)*(F_aim-math_I3))) + ! err_p = wgt*sqrt(err_p/sum((F_aim-math_I3)*(F_aim-math_I3))) - write(6,'(a,es14.7,es14.7)') 'error stress = ',err_stress/err_stress_tol - write(6,*) ' ' - write(6,'(a,es14.7)') 'max abs err F', err_f - write(6,'(a,es14.7)') 'max abs err P', err_p + ! write(6,'(a,es14.7,es14.7)') 'error stress = ',err_stress/err_stress_tol + ! write(6,*) ' ' + ! write(6,'(a,es14.7)') 'max abs err F', err_f + ! write(6,'(a,es14.7)') 'max abs err P', err_p - do k=zs,ze; do j=ys,ye; do i=xs,xe - temp = math_mul3333xx33(S_inc0,P(i,j,k,1:3,1:3)) + math_I3 - F_lambda(i,j,k,1:3,1:3) & - + F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) - f_scal(0,i,j,k) = temp(1,1) - f_scal(1,i,j,k) = temp(1,2) - f_scal(2,i,j,k) = temp(1,3) - f_scal(3,i,j,k) = temp(2,1) - f_scal(4,i,j,k) = temp(2,2) - f_scal(5,i,j,k) = temp(2,3) - f_scal(6,i,j,k) = temp(3,1) - f_scal(7,i,j,k) = temp(3,2) - f_scal(8,i,j,k) = temp(3,3) - f_scal(9,i,j,k) = F(i,j,k,1,1) - field_real(i,j,k,1,1) - f_scal(10,i,j,k) = F(i,j,k,1,2) - field_real(i,j,k,1,2) - f_scal(11,i,j,k) = F(i,j,k,1,3) - field_real(i,j,k,1,3) - f_scal(12,i,j,k) = F(i,j,k,2,1) - field_real(i,j,k,2,1) - f_scal(13,i,j,k) = F(i,j,k,2,2) - field_real(i,j,k,2,2) - f_scal(14,i,j,k) = F(i,j,k,2,3) - field_real(i,j,k,2,3) - f_scal(15,i,j,k) = F(i,j,k,3,1) - field_real(i,j,k,3,1) - f_scal(16,i,j,k) = F(i,j,k,3,2) - field_real(i,j,k,3,2) - f_scal(17,i,j,k) = F(i,j,k,3,3) - field_real(i,j,k,3,3) - enddo; enddo; enddo + ! do k=zs,ze; do j=ys,ye; do i=xs,xe + ! temp = math_mul3333xx33(S_inc0,P(i,j,k,1:3,1:3)) + math_I3 - F_lambda(i,j,k,1:3,1:3) & + ! + F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) + ! f_scal(0,i,j,k) = temp(1,1) + ! f_scal(1,i,j,k) = temp(1,2) + ! f_scal(2,i,j,k) = temp(1,3) + ! f_scal(3,i,j,k) = temp(2,1) + ! f_scal(4,i,j,k) = temp(2,2) + ! f_scal(5,i,j,k) = temp(2,3) + ! f_scal(6,i,j,k) = temp(3,1) + ! f_scal(7,i,j,k) = temp(3,2) + ! f_scal(8,i,j,k) = temp(3,3) + ! f_scal(9,i,j,k) = F(i,j,k,1,1) - field_real(i,j,k,1,1) + ! f_scal(10,i,j,k) = F(i,j,k,1,2) - field_real(i,j,k,1,2) + ! f_scal(11,i,j,k) = F(i,j,k,1,3) - field_real(i,j,k,1,3) + ! f_scal(12,i,j,k) = F(i,j,k,2,1) - field_real(i,j,k,2,1) + ! f_scal(13,i,j,k) = F(i,j,k,2,2) - field_real(i,j,k,2,2) + ! f_scal(14,i,j,k) = F(i,j,k,2,3) - field_real(i,j,k,2,3) + ! f_scal(15,i,j,k) = F(i,j,k,3,1) - field_real(i,j,k,3,1) + ! f_scal(16,i,j,k) = F(i,j,k,3,2) - field_real(i,j,k,3,2) + ! f_scal(17,i,j,k) = F(i,j,k,3,3) - field_real(i,j,k,3,3) + ! enddo; enddo; enddo - return -end subroutine FormFunctionLocal + ! return +! end subroutine FormFunctionLocal -! --------------------------------------------------------------------- -! User defined convergence check -! -subroutine SNESConverged_Interactive(snes,it,xnorm,snorm,fnorm,reason,dummy,ierr_psc) +! ! --------------------------------------------------------------------- +! ! User defined convergence check +! ! +! subroutine SNESConverged_Interactive(snes,it,xnorm,snorm,fnorm,reason,dummy,ierr_psc) - implicit none -#include + ! implicit none +! #include -! Input/output variables: - SNES snes - PetscInt it - PetscReal xnorm, snorm, fnorm - SNESConvergedReason reason - PetscObject dummy - PetscErrorCode ierr_psc +! ! Input/output variables: + ! SNES snes + ! PetscInt it + ! PetscReal xnorm, snorm, fnorm + ! SNESConvergedReason reason + ! PetscObject dummy + ! PetscErrorCode ierr_psc - err_crit = max(err_stress/err_stress_tol, & - err_f/1e-6, err_p/1e-5) - !fnorm*wgt/sqrt(sum((F_star_av-math_I3)*(F_star_av-math_I3)))/err_div_tol) + ! err_crit = max(err_stress/err_stress_tol, & + ! err_f/1e-6, err_p/1e-5) + ! !fnorm*wgt/sqrt(sum((F_star_av-math_I3)*(F_star_av-math_I3)))/err_div_tol) - if ((err_crit > 1.0_pReal .or. it < itmin) .and. it < itmax) then - reason = 0 - else - reason = 1 - endif + ! if ((err_crit > 1.0_pReal .or. it < itmin) .and. it < itmax) then + ! reason = 0 + ! else + ! reason = 1 + ! endif - return -end subroutine SNESConverged_Interactive + ! return +! end subroutine SNESConverged_Interactive -end module DAMASK_spectral_SolverAL +! end module DAMASK_spectral_SolverAL diff --git a/code/DAMASK_spectral_SolverBasic.f90 b/code/DAMASK_spectral_SolverBasic.f90 index 363b0c0b3..96ed3e9a6 100644 --- a/code/DAMASK_spectral_SolverBasic.f90 +++ b/code/DAMASK_spectral_SolverBasic.f90 @@ -1,5 +1,5 @@ !-------------------------------------------------------------------------------------------------- -!* $Id$ +! $Id$ !-------------------------------------------------------------------------------------------------- !> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH @@ -10,170 +10,194 @@ module DAMASK_spectral_SolverBasic use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran >4.6 at the moment) - use DAMASK_spectral_Utilities + use prec, only: & + pInt, & + pReal use math, only: & - math_I3, & - math_mul33x33, - - use mesh, only : & - mesh_spectral_getResolution, & - mesh_spectral_getDimension, & - math_rotate_backward33, & - math_transpose33,& - math_mul3333xx33, & - math_eigenvalues33 + math_I3 + + use DAMASK_spectral_Utilities, only: & + solutionState implicit none - character (len=*), parameter, public :: & DAMASK_spectral_SolverBasic_label = 'basic' !-------------------------------------------------------------------------------------------------- -! common pointwise data - real(pReal), dimension(:,:,:,:,:), allocatable :: F, F_lastInc, P - real(pReal), dimension(:,:,:,:), allocatable :: coordinates - real(pReal), dimension(:,:,:), allocatable :: temperature +! pointwise data + real(pReal), private, dimension(:,:,:,:,:), allocatable :: F, F_lastInc, P + real(pReal), private, dimension(:,:,:,:), allocatable :: coordinates + real(pReal), private, dimension(:,:,:), allocatable :: temperature !-------------------------------------------------------------------------------------------------- ! stress, stiffness and compliance average etc. - real(pReal), dimension(3,3) :: & + real(pReal), private, dimension(3,3) :: & F_aim = math_I3, & F_aim_lastInc = math_I3 - real(pReal), dimension(3,3,3,3) :: & + real(pReal), private,dimension(3,3,3,3) :: & C = 0.0_pReal - contains - subroutine basic_init() +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields and fills them with data, potentially from restart info +!-------------------------------------------------------------------------------------------------- +subroutine basic_init() - 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 - implicit none - integer(pInt) :: i,j,k + use DAMASK_spectral_Utilities, only: & + Utilities_init, & + Utilities_constitutiveResponse, & + Utilities_updateGamma, & + debugrestart + + use mesh, only: & + res, & + geomdim - real(pReal), dimension(3,3) :: temp33_Real + implicit none + + integer(pInt) :: i,j,k + real(pReal), dimension(3,3) :: temp33_Real - write(6,'(a)') '' - write(6,'(a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' - write(6,'(a)') ' $Id$' + + call Utilities_Init() + write(6,'(a)') '' + write(6,'(a)') ' <<<+- DAMASK_spectral_solverBasic init -+>>>' + write(6,'(a)') ' $Id$' #include "compilation_info.f90" - write(6,'(a)') '' - call Utilities_Init() + write(6,'(a)') '' + - allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) - allocate (P ( res(1), res(2),res(3),3,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) + allocate (F ( res(1), res(2),res(3),3,3), source = 0.0_pReal) + allocate (F_lastInc ( res(1), res(2),res(3),3,3), source = 0.0_pReal) + allocate (P ( res(1), res(2),res(3),3,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) !-------------------------------------------------------------------------------------------------- ! init fields - if (restartInc == 1_pInt) then ! no deformation (no restart) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = math_I3 - F_lastInc(i,j,k,1:3,1:3) = math_I3 - coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & - - geomdim/real(2_pInt*res,pReal) - enddo; enddo; enddo - elseif (restartInc > 1_pInt) then ! using old values from file - if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& + if (restartInc == 1_pInt) then ! no deformation (no restart) + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + F(i,j,k,1:3,1:3) = math_I3 + F_lastInc(i,j,k,1:3,1:3) = math_I3 + coordinates(i,j,k,1:3) = geomdim/real(res,pReal)*real([i,j,k],pReal) & + - geomdim/real(2_pInt*res,pReal) + enddo; enddo; enddo + elseif (restartInc > 1_pInt) then ! using old values from file + if (debugRestart) write(6,'(a,i6,a)') 'Reading values of increment ',& restartInc - 1_pInt,' from file' - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad',& trim(getSolverJobName()),size(F)) - read (777,rec=1) F - close (777) - call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& + read (777,rec=1) F + close (777) + call IO_read_jobBinaryFile(777,'convergedSpectralDefgrad_lastInc',& trim(getSolverJobName()),size(F_lastInc)) - read (777,rec=1) F_lastInc - close (777) - call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) - read (777,rec=1) F_aim - close (777) - call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) - read (777,rec=1) F_aim_lastInc - close (777) + read (777,rec=1) F_lastInc + close (777) + call IO_read_jobBinaryFile(777,'F_aim',trim(getSolverJobName()),size(F_aim)) + read (777,rec=1) F_aim + close (777) + call IO_read_jobBinaryFile(777,'F_aim_lastInc',trim(getSolverJobName()),size(F_aim_lastInc)) + read (777,rec=1) F_aim_lastInc + close (777) - coordinates = 0.0 ! change it later!!! - endif + coordinates = 0.0 ! change it later!!! + endif - call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& - P,C,temp33_Real,.false.,math_I3) + call Utilities_constitutiveResponse(coordinates,F,F_lastInc,temperature,0.0_pReal,& + P,C,temp33_Real,.false.,math_I3) !-------------------------------------------------------------------------------------------------- ! reference stiffness - if (restartInc == 1_pInt) then - call IO_write_jobBinaryFile(777,'C_ref',size(C)) - write (777,rec=1) C - close(777) - elseif (restartInc > 1_pInt) then - call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) - read (777,rec=1) C - close (777) - endif + if (restartInc == 1_pInt) then + call IO_write_jobBinaryFile(777,'C_ref',size(C)) + write (777,rec=1) C + close(777) + elseif (restartInc > 1_pInt) then + call IO_read_jobBinaryFile(777,'C_ref',trim(getSolverJobName()),size(C)) + read (777,rec=1) C + close (777) + endif - call Utilities_updateGamma(C) + call Utilities_updateGamma(C) - end subroutine basic_init - -type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F_BC,mask_stressVector,velgrad,rotation_BC) +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) use numerics, only: & itmax, & itmin, & update_gamma - + use math, only: & + math_mul33x33 ,& + math_rotate_backward33, & + math_transpose33, & + math_mul3333xx33, & + deformed_fft + use mesh, only: & + res,& + geomdim use IO, only: & IO_write_JobBinaryFile + use DAMASK_spectral_Utilities, only: & + boundaryCondition, & + field_real, & + Utilities_forwardField, & + Utilities_maskedCompliance, & + Utilities_forwardFFT, & + Utilities_divergenceRMS, & + Utilities_fourierConvolution, & + Utilities_backwardFFT, & + Utilities_updateGamma, & + Utilities_constitutiveResponse + use FEsolving, only: & restartWrite implicit none - !-------------------------------------------------------------------------------------------------- ! input data for solution - - real(pReal), intent(in) :: timeinc, timeinc_old - real(pReal), intent(in) :: guessmode - logical, intent(in) :: velgrad - real(pReal), dimension(3,3), intent(in) :: P_BC,F_BC,rotation_BC - logical, dimension(9), intent(in) :: mask_stressVector + real(pReal), intent(in) :: timeinc, timeinc_old, temperature_bc, guessmode + type(boundaryCondition), intent(in) :: P_BC,F_BC + real(pReal), dimension(3,3), intent(in) :: rotation_BC -!-------------------------------------------------------------------------------------------------- -! loop variables, convergence etc. - real(pReal), dimension(3,3), parameter :: ones = 1.0_pReal, zeroes = 0.0_pReal - real(pReal), dimension(3,3) :: temp33_Real + real(pReal), dimension(3,3,3,3) :: S - real(pReal), dimension(3,3) :: mask_stress, & - mask_defgrad, & - deltaF_aim, & + real(pReal), dimension(3,3) :: deltaF_aim, & F_aim_lab, & F_aim_lab_lastIter, & P_av +!-------------------------------------------------------------------------------------------------- +! loop variables, convergence etc. real(pReal) :: err_div, err_stress - integer(pInt) :: iter - integer(pInt) :: i, j, k + integer(pInt) :: iter, row, column, i, j, k logical :: ForwardData - real(pReal) :: defgradDet - real(pReal) :: defgradDetMax, defgradDetMin + real(pReal) :: defgradDet, defgradDetMax, defgradDetMin + real(pReal), dimension(3,3) :: temp33_Real - mask_stress = merge(ones,zeroes,reshape(mask_stressVector,[3,3])) - mask_defgrad = merge(zeroes,ones,reshape(mask_stressVector,[3,3])) - +!-------------------------------------------------------------------------------------------------- +! restart information for spectral solver if (restartWrite) then write(6,'(a)') 'writing converged results for restart' - call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) ! writing deformation gradient field to file + call IO_write_jobBinaryFile(777,'convergedSpectralDefgrad',size(F_lastInc)) write (777,rec=1) F_LastInc close (777) call IO_write_jobBinaryFile(777,'C',size(C)) @@ -181,18 +205,16 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F close(777) endif - ForwardData = .True. - if (velgrad) then ! calculate deltaF_aim from given L and current F - deltaF_aim = timeinc * mask_defgrad * math_mul33x33(F_BC, F_aim) - else ! deltaF_aim = fDot *timeinc where applicable - deltaF_aim = timeinc * mask_defgrad * F_BC - endif - !-------------------------------------------------------------------------------------------------- ! 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 + endif temp33_Real = F_aim F_aim = F_aim & - + guessmode * mask_stress * (F_aim - F_aim_lastInc)*timeinc/timeinc_old & + + 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 @@ -200,21 +222,17 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F !-------------------------------------------------------------------------------------------------- ! update local deformation gradient and coordinates deltaF_aim = math_rotate_backward33(deltaF_aim,rotation_BC) - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - temp33_Real = F(i,j,k,1:3,1:3) - F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon - + guessmode * (F(i,j,k,1:3,1:3) - F_lastInc(i,j,k,1:3,1:3))*timeinc/timeinc_old& ! guessing... - + (1.0_pReal-guessmode) * deltaF_aim ! if not guessing, use prescribed average deformation where applicable - F_lastInc(i,j,k,1:3,1:3) = temp33_Real - enddo; enddo; enddo - call deformed_fft(res,geomdim,math_rotate_backward33(F_aim,rotation_BC),& ! calculate current coordinates - 1.0_pReal,F_lastInc,coordinates) + 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) - iter = 0_pInt - S = Utilities_maskedCompliance(rotation_BC,mask_stressVector,C) +!-------------------------------------------------------------------------------------------------- +! update stiffness (and gamma operator) + S = Utilities_maskedCompliance(rotation_BC,P_BC%maskLogical,C) if (update_gamma) call Utilities_updateGamma(C) - - convergenceLoop: do + + iter = 0_pInt + ForwardData = .True. + convergenceLoop: do while(iter < itmax) iter = iter + 1_pInt !-------------------------------------------------------------------------------------------------- @@ -222,78 +240,94 @@ type(solutionState) function basic_solution(guessmode,timeinc,timeinc_old,P_BC,F write(6,'(a)') '' write(6,'(a)') '==================================================================' write(6,'(3(a,i6.6))') ' Iter. ',itmin,' < ',iter,' < ',itmax + 1_pInt - write(6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'deformation gradient aim =',& - math_transpose33(F_aim) + 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 -print*, 'FLast 111', F_lastInc(1,1,1,1:3,1:3) call Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,timeinc,& - P,C,P_av,ForwardData,rotation_BC) + P,C,P_av,ForwardData,rotation_BC) ForwardData = .False. !-------------------------------------------------------------------------------------------------- ! stress BC handling - F_aim = F_aim - math_mul3333xx33(S, ((P_av - P_BC))) !S = 0.0 for no bc - err_stress = maxval(mask_stress * (P_av - P_BC)) ! 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(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 - field_real = 0.0_pReal - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3) = P - call Utilities_forwardFFT() - err_div = Utilities_divergenceRMS() - call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) - call Utilities_backwardFFT() +! 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) = P + call Utilities_forwardFFT() + err_div = Utilities_divergenceRMS() + call Utilities_fourierConvolution(F_aim_lab_lastIter - F_aim_lab) + call Utilities_backwardFFT() + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization + enddo; enddo; enddo - temp33_real =0.0_pReal - - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) - F(i,j,k,1:3,1:3) = F(i,j,k,1:3,1:3) - field_real(i,j,k,1:3,1:3) ! F(x)^(n+1) = F(x)^(n) + correction; *wgt: correcting for missing normalization - temp33_real = temp33_real + field_real(i,j,k,1:3,1:3) - enddo; enddo; enddo - if (Convergenced(err_div,P_av,err_stress,P_av,iter)) exit - + basic_solution%converged = basic_Convergeced(err_div,P_av,err_stress,P_av) + + if (basic_solution%converged .and. iter > itmin) exit enddo convergenceLoop end function basic_solution -logical function Convergenced(err_div,P_av,err_stress,P_av2,iter) + +!-------------------------------------------------------------------------------------------------- +!> @brief convergence check for basic scheme based on div of P and deviation from stress aim +!-------------------------------------------------------------------------------------------------- +logical function basic_Convergeced(err_div,pAvgDiv,err_stress,pAvgStress) use numerics, only: & - itmax, & - itmin, & - err_div_tol, & - err_stress_tolrel, & - err_stress_tolabs + itmin, & + err_div_tol, & + err_stress_tolrel, & + err_stress_tolabs + + use math, only: & + math_mul33x33, & + math_eigenvalues33, & + math_transpose33 implicit none - real(pReal), dimension(3,3) :: P_av, P_av2 - real(pReal) :: err_div, err_stress, field_av_L2 - integer(pInt) :: iter + real(pReal), dimension(3,3), intent(in) :: & + pAvgDiv,& + pAvgStress + + real(pReal), intent(in) :: & + err_div, & + err_stress + + real(pReal) :: & + err_stress_tol, & + pAvgDivL2 + + - field_av_L2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(P_av,& ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) - math_transpose33(P_av))))) - Convergenced = (iter < itmax) .and. (iter > itmin) .and. & - all([err_div/field_av_L2/err_div_tol,& - err_stress/min(maxval(abs(P_av2))*err_stress_tolrel,err_stress_tolabs)] < 1.0_pReal) + pAvgDivL2 = sqrt(maxval(math_eigenvalues33(math_mul33x33(pAvgDiv,math_transpose33(pAvgDiv))))) ! L_2 norm of average stress (http://mathworld.wolfram.com/SpectralNorm.html) + err_stress_tol = min(maxval(abs(pAvgStress))*err_stress_tolrel,err_stress_tolabs) + + basic_Convergeced = all([ err_div/pAvgDivL2/err_div_tol,& + err_stress/err_stress_tol ] < 1.0_pReal) - - write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/min(maxval(abs(P_av2))*err_stress_tolrel,err_stress_tolabs), & - ' (',err_stress,' Pa)' - write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/field_av_L2/err_div_tol,& - ' (',err_div,' N/m³)' -end function Convergenced + write(6,'(a,f6.2,a,es11.4,a)') 'error divergence = ', err_div/pAvgDivL2/err_div_tol,& + ' (',err_div,' N/m³)' + write(6,'(a,f6.2,a,es11.4,a)') 'error stress = ', err_stress/err_stress_tol, & + ' (',err_stress,' Pa)' + +end function basic_Convergeced subroutine basic_destroy() - -implicit none -call Utilities_destroy() + + use DAMASK_spectral_Utilities, only: & + Utilities_destroy + + implicit none + call Utilities_destroy() end subroutine basic_destroy diff --git a/code/DAMASK_spectral_Utilities.f90 b/code/DAMASK_spectral_Utilities.f90 index 8987b96b3..3cbabb697 100644 --- a/code/DAMASK_spectral_Utilities.f90 +++ b/code/DAMASK_spectral_Utilities.f90 @@ -11,7 +11,13 @@ module DAMASK_spectral_Utilities use prec, only: & pReal, & pInt - + use mesh, only : & + res, & + res1_red, & + geomdim, & + mesh_NcpElems, & + wgt + use math use IO, only: & @@ -21,16 +27,17 @@ module DAMASK_spectral_Utilities !-------------------------------------------------------------------------------------------------- ! variables storing information for spectral method and FFTW - type(C_PTR) , private :: plan_forward, plan_backward ! plans for fftw - real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method - real(pReal), private, dimension(3,3,3,3) :: C_ref - real(pReal), private, dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator - real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real - complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier + type(C_PTR), private :: plan_forward, plan_backward ! plans for fftw + real(pReal), private, dimension(:,:,:,:,:,:,:), allocatable :: gamma_hat ! gamma operator (field) for spectral method + real(pReal), private, dimension(:,:,:,:), allocatable :: xi ! wave vector field for divergence and for gamma operator + complex(pReal),private, dimension(:,:,:,:,:), pointer :: field_fourier + real(pReal), private, dimension(3,3,3,3) :: C_ref + + real(pReal), public, dimension(:,:,:,:,:), pointer :: field_real !-------------------------------------------------------------------------------------------------- ! debug fftw - type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back + type(C_PTR), private :: plan_scalarField_forth, plan_scalarField_back complex(pReal),private, dimension(:,:,:), pointer :: scalarField_real complex(pReal),private, dimension(:,:,:), pointer :: scalarField_fourier @@ -41,40 +48,33 @@ module DAMASK_spectral_Utilities complex(pReal), private, dimension(:,:,:,:), pointer :: divergence_fourier real(pReal), dimension(:,:,:,:), allocatable :: divergence_post - type BC_type - real(pReal), dimension(3,3) :: values - real(pReal), dimension(3,3) :: maskFloat - logical, dimension(3,3) :: maskLogical - character(20) :: myType - end type BC_type - !-------------------------------------------------------------------------------------------------- !variables controlling debugging - logical :: debugGeneral, debugDivergence, debugRestart, debugFFTW - - - real(pReal), dimension(3) :: geomdim = 0.0_pReal, virt_dim = 0.0_pReal ! physical dimension of volume element per direction - integer(pInt), dimension(3) :: res = 1_pInt - real(pReal) :: wgt - integer(pInt) :: res1_red, Npoints + logical,public :: debugGeneral, debugDivergence, debugRestart, debugFFTW !-------------------------------------------------------------------------------------------------- -! solution state +! derived types type solutionState logical :: converged = .false. logical :: regrid = .false. logical :: term_ill = .false. end type solutionState + + type boundaryCondition + 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 + contains +!-------------------------------------------------------------------------------------------------- +!> @brief allocates all neccessary fields, sets debug flags, create plans for fftw +!-------------------------------------------------------------------------------------------------- subroutine Utilities_init() - - use mesh, only : & - mesh_spectral_getResolution, & - mesh_spectral_getDimension - - use numerics, only: & - divergence_correction, & + + use numerics, only: & DAMASK_NumThreadsInt, & fftw_planner_flag, & fftw_timelimit, & @@ -88,11 +88,13 @@ subroutine Utilities_init() debug_spectralRestart, & debug_spectralFFTW + use mesh, only : & + virt_dim + implicit none - - integer(pInt) :: i, j, k, ierr + integer(pInt) :: i, j, k integer(pInt), dimension(3) :: k_s - + !$ integer(pInt) :: ierr type(C_PTR) :: tensorField ! field in real and fourier space type(C_PTR) :: scalarField_realC, scalarField_fourierC type(C_PTR) :: divergence @@ -111,15 +113,8 @@ subroutine Utilities_init() debugRestart = iand(debug_level(debug_spectral),debug_spectralRestart) /= 0 debugFFTW = iand(debug_level(debug_spectral),debug_spectralFFTW) /= 0 -!################################################################################################## -! initialization -!################################################################################################## - res = mesh_spectral_getResolution() - geomdim = mesh_spectral_getDimension() - res1_red = res(1)/2_pInt + 1_pInt - Npoints = res(1)*res(2)*res(3) - wgt = 1.0/real(Npoints,pReal) - +!-------------------------------------------------------------------------------------------------- +! allocation allocate (xi (3,res1_red,res(2),res(3)), source = 0.0_pReal) ! start out isothermally tensorField = fftw_alloc_complex(int(res1_red*res(2)*res(3)*9_pInt,C_SIZE_T)) ! allocate continous data using a C function, C_SIZE_T is of type integer(8) call c_f_pointer(tensorField, field_real, [ res(1)+2_pInt,res(2),res(3),3,3]) ! place a pointer for a real representation on tensorField @@ -178,14 +173,6 @@ subroutine Utilities_init() !-------------------------------------------------------------------------------------------------- ! calculation of discrete angular frequencies, ordered as in FFTW (wrap around) - if (divergence_correction) then - do i = 1_pInt, 3_pInt - if (i /= minloc(geomdim,1) .and. i /= maxloc(geomdim,1)) virt_dim = geomdim/geomdim(i) - enddo - else - virt_dim = geomdim - endif - do k = 1_pInt, res(3) k_s(3) = k - 1_pInt if(k > res(3)/2_pInt + 1_pInt) k_s(3) = k_s(3) - res(3) @@ -202,59 +189,65 @@ subroutine Utilities_init() else ! precalculation of gamma_hat field allocate (gamma_hat(res1_red ,res(2),res(3),3,3,3,3), source =0.0_pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 endif + end subroutine Utilities_init +!-------------------------------------------------------------------------------------------------- +!> @brief updates references stiffness and potentially precalculated gamma operator +!-------------------------------------------------------------------------------------------------- subroutine Utilities_updateGamma(C) - use numerics, only: & - memory_efficient + use numerics, only: & + memory_efficient - implicit none + implicit none - real(pReal), dimension(3,3,3,3) :: C - real(pReal), dimension(3,3) :: temp33_Real, xiDyad - integer(pInt) :: i, j, k, l, m, n, o + real(pReal), dimension(3,3,3,3), intent(in) :: C + real(pReal), dimension(3,3) :: temp33_Real, xiDyad + integer(pInt) :: i, j, k, l, m, n, o - C_ref = C - if(.not. memory_efficient) then - do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red - if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) - forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & - temp33_Real(l,m) = sum(C_ref(l,m,1:3,1:3)*xiDyad) - temp33_Real = math_inv33(temp33_Real) - forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& - gamma_hat(i,j,k, l,m,n,o) = temp33_Real(l,n)*xiDyad(m,o) - endif - enddo; enddo; enddo - gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 - endif - + C_ref = C + if(.not. memory_efficient) then + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red + if(any([i,j,k] /= 1_pInt)) then ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + xiDyad(l,m) = xi(l, i,j,k)*xi(m, i,j,k) + forall(l = 1_pInt:3_pInt, m = 1_pInt:3_pInt) & + temp33_Real(l,m) = sum(C_ref(l,m,1:3,1:3)*xiDyad) + temp33_Real = math_inv33(temp33_Real) + forall(l=1_pInt:3_pInt, m=1_pInt:3_pInt, n=1_pInt:3_pInt, o=1_pInt:3_pInt)& + gamma_hat(i,j,k, l,m,n,o) = temp33_Real(l,n)*xiDyad(m,o) + endif + enddo; enddo; enddo + gamma_hat(1,1,1, 1:3,1:3,1:3,1:3) = 0.0_pReal ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + endif end subroutine Utilities_updateGamma -subroutine Utilities_forwardFFT() - - implicit none +!-------------------------------------------------------------------------------------------------- +!> @brief forward FFT of data in field_real to field_fourier with highest freqs. removed +!-------------------------------------------------------------------------------------------------- +subroutine Utilities_forwardFFT(row,column) + use mesh, only : & + virt_dim - integer(pInt) :: row, column + implicit none + integer(pInt), intent(in), optional :: row, column !-------------------------------------------------------------------------------------------------- ! copy one component of the stress field to to a single FT and check for mismatch if (debugFFTW) then + if (.not. present(row) .or. .not. present(column)) stop scalarField_real(1:res(1),1:res(2),1:res(3)) =& ! store the selected component cmplx(field_real(1:res(1),1:res(2),1:res(3),row,column),0.0_pReal,pReal) - endif !-------------------------------------------------------------------------------------------------- ! call function to calculate divergence from math (for post processing) to check results if (debugDivergence) & - call divergence_fft(res,virt_dim,3_pInt,& - field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),divergence_post) + call divergence_fft(res,virt_dim,3_pInt,field_real(1:res(1),1:res(2),1:res(3),1:3,1:3),divergence_post) !-------------------------------------------------------------------------------------------------- -! doing the FT because it simplifies calculation of average stress in real space also +! doing the FT call fftw_execute_dft_r2c(plan_forward,field_real,field_fourier) !-------------------------------------------------------------------------------------------------- @@ -281,17 +274,16 @@ subroutine Utilities_forwardFFT() = cmplx(0.0_pReal,0.0_pReal,pReal) end subroutine Utilities_forwardFFT -subroutine Utilities_backwardFFT() +subroutine Utilities_backwardFFT(row,column) implicit none - integer(pInt) :: row, column, i, j, k, m, n + integer(pInt), intent(in), optional :: row, column + integer(pInt) :: i, j, k, m, n !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 inverse FT results if (debugFFTW) then - row = 3 ! (mod(totalIncsCounter+iter-2_pInt,9_pInt))/3_pInt + 1_pInt ! go through the elements of the tensors, controlled by totalIncsCounter and iter, starting at 1 - column = 3 !(mod(totalIncsCounter+iter-2_pInt,3_pInt)) + 1_pInt do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res1_red scalarField_fourier(i,j,k) = field_fourier(i,j,k,row,column) enddo; enddo; enddo @@ -308,11 +300,11 @@ subroutine Utilities_backwardFFT() m = m -1_pInt enddo; enddo endif + !-------------------------------------------------------------------------------------------------- ! doing the inverse FT - print*, 'field fourier 111', field_fourier(1,1,1,1:3,1:3) call fftw_execute_dft_c2r(plan_backward,field_fourier,field_real) ! back transform of fluct deformation gradient - print*, 'field real 111', field_real(1,1,1,1:3,1:3) + !-------------------------------------------------------------------------------------------------- ! comparing 1 and 3x3 inverse FT results if (debugFFTW) then @@ -324,6 +316,7 @@ subroutine Utilities_backwardFFT() real(scalarField_real(1:res(1),1:res(2),1:res(3)))) endif field_real = field_real * wgt + end subroutine Utilities_backwardFFT @@ -370,7 +363,7 @@ subroutine Utilities_fourierConvolution(fieldAim) enddo; enddo; enddo endif - field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(Npoints,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 + field_fourier(1,1,1,1:3,1:3) = cmplx(fieldAim*real(mesh_NcpElems,pReal),0.0_pReal,pReal) ! singular point at xi=(0.0,0.0,0.0) i.e. i=j=k=1 end subroutine Utilities_fourierConvolution @@ -390,7 +383,6 @@ real(pReal) function Utilities_divergenceRMS() !-------------------------------------------------------------------------------------------------- ! calculating RMS divergence criterion in Fourier space - Utilities_divergenceRMS = 0.0_pReal do k = 1_pInt, res(3); do j = 1_pInt, res(2) do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice. @@ -440,11 +432,11 @@ real(pReal) function Utilities_divergenceRMS() err_real_div_max = sqrt( err_real_div_max) ! max in real space err_div_max = sqrt( err_div_max) ! max in Fourier space - write(6,'(a,es11.4)') 'error divergence FT RMS = ',err_div_RMS - write(6,'(a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS - write(6,'(a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS - write(6,'(a,es11.4)') 'error divergence FT max = ',err_div_max - write(6,'(a,es11.4)') 'error divergence Real max = ',err_real_div_max + write(6,'(1x,a,es11.4)') 'error divergence FT RMS = ',err_div_RMS + write(6,'(1x,a,es11.4)') 'error divergence Real RMS = ',err_real_div_RMS + write(6,'(1x,a,es11.4)') 'error divergence post RMS = ',err_post_div_RMS + write(6,'(1x,a,es11.4)') 'error divergence FT max = ',err_div_max + write(6,'(1x,a,es11.4)') 'error divergence Real max = ',err_real_div_max endif end function Utilities_divergenceRMS @@ -536,8 +528,6 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti CPFEM_mode = 2_pInt endif - write(6,'(a)') '' - write(6,'(a)') '... update stress field P(F) .....................................' 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 @@ -564,13 +554,30 @@ subroutine Utilities_constitutiveResponse(coordinates,F_lastInc,F,temperature,ti restartWrite = .false. P_av_lab = P_av_lab * wgt P_av = math_rotate_forward33(P_av_lab,rotation_BC) - write (6,'(a,/,3(3(f12.7,1x)/))',advance='no') 'Piola-Kirchhoff stress / MPa =',& - math_transpose33(P_av)/1.e6_pReal + 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) + + real(pReal), intent(in), dimension(3,3) :: delta_aim + + real(pReal), intent(in) :: timeinc, timeinc_old, guessmode + real(pReal), intent(inout), dimension(res(1),res(2),res(3),3,3) :: field_lastInc,field + integer(pInt) :: i,j,k + real(pReal), dimension(3,3) :: temp33_real + + do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1) + temp33_Real = Field(i,j,k,1:3,1:3) + Field(i,j,k,1:3,1:3) = Field(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon + + guessmode * (field(i,j,k,1:3,1:3) - Field_lastInc(i,j,k,1:3,1:3))*timeinc/timeinc_old& ! guessing... + + (1.0_pReal-guessmode) * delta_aim ! if not guessing, use prescribed average deformation where applicable + Field_lastInc(i,j,k,1:3,1:3) = temp33_Real + enddo; enddo; enddo + end subroutine Utilities_forwardField + subroutine Utilities_destroy() implicit none diff --git a/code/DAMASK_spectral_interface.f90 b/code/DAMASK_spectral_interface.f90 index 4a77b623e..302a89db2 100644 --- a/code/DAMASK_spectral_interface.f90 +++ b/code/DAMASK_spectral_interface.f90 @@ -20,6 +20,7 @@ !* $Id$ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief Interfacing between the spectral solver and the material subroutines provided !! by DAMASK !-------------------------------------------------------------------------------------------------- @@ -100,48 +101,48 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) tag = IO_lc(IO_stringValue(commandLine,positions,i)) ! extract key select case(tag) case ('-h','--help') - write(6,'(a)') '#############################################################' - write(6,'(a)') 'DAMASK spectral:' - write(6,'(a)') 'The spectral method boundary value problem solver for' - write(6,'(a)') 'the Duesseldorf Advanced Material Simulation Kit' - write(6,'(a)') '#############################################################' - write(6,'(a)') 'Valid command line switches:' - write(6,'(a)') ' --geom (-g, --geometry)' - write(6,'(a)') ' --load (-l, --loadcase)' - write(6,'(a)') ' --restart (-r, --rs)' - write(6,'(a)') ' --regrid (--rg)' - write(6,'(a)') ' --help (-h)' - write(6,'(a)') ' ' - write(6,'(a)') 'Mandatory Arguments:' - write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load' - write(6,'(a)') ' "PathToLoadFile" will be the working directory.' - write(6,'(a)') ' Make sure the file "material.config" exists in the working' - write(6,'(a)') ' directory' - write(6,'(a)') ' For further configuration place "numerics.config"' - write(6,'(a)') ' and "numerics.config" in that directory.' - write(6,'(a)') ' ' - write(6,'(a)') ' --geom PathToGeomFile/NameOfGeom.geom' - write(6,'(a)') ' ' - write(6,'(a)') 'Optional Argument:' - write(6,'(a)') ' --restart XX' - write(6,'(a)') ' Reads in total increment No. XX-1 and continous to' - write(6,'(a)') ' calculate total increment No. XX.' - write(6,'(a)') ' Appends to existing results file ' - write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' - write(6,'(a)') ' Works only if the restart information for total increment' - write(6,'(a)') ' No. XX-1 is available in the working directory.' - write(6,'(a)') ' ' - write(6,'(a)') ' --regrid XX' - write(6,'(a)') ' Reads in total increment No. XX-1 and continous to' - write(6,'(a)') ' calculate total increment No. XX.' - write(6,'(a)') ' Attention: Overwrites existing results file ' - write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' - write(6,'(a)') ' Works only if the restart information for total increment' - write(6,'(a)') ' No. XX-1 is available in the working directory.' - write(6,'(a)') 'Help:' - write(6,'(a)') ' --help' - write(6,'(a)') ' Prints this message and exits' - write(6,'(a)') ' ' + write(6,'(a)') ' #############################################################' + write(6,'(a)') ' DAMASK spectral:' + write(6,'(a)') ' The spectral method boundary value problem solver for' + write(6,'(a)') ' the Duesseldorf Advanced Material Simulation Kit' + write(6,'(a)') ' #############################################################' + write(6,'(a)') ' Valid command line switches:' + write(6,'(a)') ' --geom (-g, --geometry)' + write(6,'(a)') ' --load (-l, --loadcase)' + write(6,'(a)') ' --restart (-r, --rs)' + write(6,'(a)') ' --regrid (--rg)' + write(6,'(a)') ' --help (-h)' + write(6,'(a)') ' ' + write(6,'(a)') ' Mandatory Arguments:' + write(6,'(a)') ' --load PathToLoadFile/NameOfLoadFile.load' + write(6,'(a)') ' "PathToLoadFile" will be the working directory.' + write(6,'(a)') ' Make sure the file "material.config" exists in the working' + write(6,'(a)') ' directory' + write(6,'(a)') ' For further configuration place "numerics.config"' + write(6,'(a)') ' and "numerics.config" in that directory.' + write(6,'(a)') ' ' + write(6,'(a)') ' --geom PathToGeomFile/NameOfGeom.geom' + write(6,'(a)') ' ' + write(6,'(a)') ' Optional Argument:' + write(6,'(a)') ' --restart XX' + write(6,'(a)') ' Reads in total increment No. XX-1 and continous to' + write(6,'(a)') ' calculate total increment No. XX.' + write(6,'(a)') ' Appends to existing results file ' + write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' + write(6,'(a)') ' Works only if the restart information for total increment' + write(6,'(a)') ' No. XX-1 is available in the working directory.' + write(6,'(a)') ' ' + write(6,'(a)') ' --regrid XX' + write(6,'(a)') ' Reads in total increment No. XX-1 and continous to' + write(6,'(a)') ' calculate total increment No. XX.' + write(6,'(a)') ' Attention: Overwrites existing results file ' + write(6,'(a)') ' "NameOfGeom_NameOfLoadFile.spectralOut".' + write(6,'(a)') ' Works only if the restart information for total increment' + write(6,'(a)') ' No. XX-1 is available in the working directory.' + write(6,'(a)') ' Help:' + write(6,'(a)') ' --help' + write(6,'(a)') ' Prints this message and exits' + write(6,'(a)') ' ' call quit(0_pInt) ! normal Termination case ('-l', '--load', '--loadcase') loadcaseParameter = IO_stringValue(commandLine,positions,i+1_pInt) @@ -160,7 +161,7 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) endif if (.not. (gotLoadCase .and. gotGeometry)) then - write(6,'(a)') 'Please specify Geometry AND Load Case' + write(6,'(a)') ' Please specify Geometry AND Load Case' call quit(1_pInt) endif @@ -171,21 +172,21 @@ subroutine DAMASK_interface_init(loadCaseParameterIn,geometryParameterIn) call get_environment_variable('USER',userName) call date_and_time(values = dateAndTime) - write(6,'(a,2(i2.2,a),i4.4)') 'Date: ',dateAndTime(3),'/',& - dateAndTime(2),'/',& - dateAndTime(1) - write(6,'(a,2(i2.2,a),i2.2)') 'Time: ',dateAndTime(5),':',& - dateAndTime(6),':',& - dateAndTime(7) - write(6,'(a,a)') 'Host name: ', trim(hostName) - write(6,'(a,a)') 'User name: ', trim(userName) - write(6,'(a,a)') 'Path separator: ', getPathSep() - write(6,'(a,a)') 'Command line call: ', trim(commandLine) - write(6,'(a,a)') 'Geometry parameter: ', trim(geometryParameter) - write(6,'(a,a)') 'Loadcase parameter: ', trim(loadcaseParameter) + write(6,'(a,2(i2.2,a),i4.4)') ' Date: ',dateAndTime(3),'/',& + dateAndTime(2),'/',& + dateAndTime(1) + write(6,'(a,2(i2.2,a),i2.2)') ' Time: ',dateAndTime(5),':',& + dateAndTime(6),':',& + dateAndTime(7) + write(6,'(a,a)') ' Host name: ', trim(hostName) + write(6,'(a,a)') ' User name: ', trim(userName) + write(6,'(a,a)') ' Path separator: ', getPathSep() + write(6,'(a,a)') ' Command line call: ', trim(commandLine) + write(6,'(a,a)') ' Geometry parameter: ', trim(geometryParameter) + write(6,'(a,a)') ' Loadcase parameter: ', trim(loadcaseParameter) if (SpectralRestartInc > 1_pInt) write(6,'(a,i6.6)') & - 'Restart at increment: ', spectralRestartInc - write(6,'(a,l1)') 'Append to result file: ', appendToOutFile + ' Restart at increment: ', spectralRestartInc + write(6,'(a,l1)') ' Append to result file: ', appendToOutFile end subroutine DAMASK_interface_init diff --git a/code/constitutive_none.f90 b/code/constitutive_none.f90 index ba35ee053..3d26c1aeb 100644 --- a/code/constitutive_none.f90 +++ b/code/constitutive_none.f90 @@ -17,7 +17,7 @@ ! along with DAMASK. If not, see . ! !############################################################## -!* $Id:$ +!* $Id$ !***************************************************** !* Module: CONSTITUTIVE_J2 * !***************************************************** @@ -114,7 +114,7 @@ subroutine constitutive_none_init(myFile) !$OMP CRITICAL (write2out) write(6,*) write(6,*) '<<<+- constitutive_',trim(constitutive_none_label),' init -+>>>' - write(6,*) '$Id:$' + write(6,*) '$Id$' #include "compilation_info.f90" !$OMP END CRITICAL (write2out)